Merge pull request #38 from danielparton/master

Ability to fix protonation states of specific residues

.travis.yml

@@ -34,14 +34,7 @@ install:
   - export CXX=g++
 
 script:
-  - conda config --add channels https://conda.binstar.org/omnia
-  - conda build devtools/conda-recipe
-  - conda install --yes --use-local ensembler-dev
-  # - source devtools/travis-ci/test.sh
-  - conda install --yes nose
-  - pushd .; cd /
-  - nosetests ensembler -v --exe -a modeller
-  - popd
+  - source devtools/travis-ci/test.sh
 
 after_success:
   - echo "after_success"

devtools/conda-recipe/meta.yaml

@@ -18,12 +18,12 @@ requirements:
 
   run:
     - python
-    - modeller
+    - modeller ==9.15
     - mdtraj
     - msmbuilder
     - biopython
     - openmm
-    - pdbfixer
+    - pdbfixer ==1.2
     # - numpy
     - lxml
     - pyyaml

devtools/travis-ci/test.sh

@@ -1,10 +1,16 @@
-# This runs unit tests
+# This installs the program and runs unit tests
+set -e
+conda config --add channels https://conda.anaconda.org/omnia
+conda build devtools/conda-recipe
+conda install --yes --use-local ensembler-dev
+conda install --yes nose
+pushd .; cd /
+nosetests ensembler -v --exe -a unit
 
 if [[ "$TRAVIS_PULL_REQUEST" == "false" ]]; then
-    conda install --yes nose
-    pushd .; cd /
     nosetests ensembler -v --exe -a modeller
-    popd
 else
     echo "This is a pull request. Secure environment variables are not available, so will not attempt to run Modeller tests."
 fi
+
+popd

ensembler/cli_commands/package_models.py

@@ -54,6 +54,7 @@
 helpstring = '\n\n'.join([helpstring_header, '\n\n'.join(helpstring_unique_options), '\n\n'.join(helpstring_nonunique_options)])
 docopt_helpstring = '\n\n'.join(helpstring_unique_options)
 
+
 def dispatch(args):
     if args['--package_for']:
         package_for = args['--package_for']

ensembler/cli_commands/refine_implicit.py

@@ -122,6 +122,6 @@ def dispatch(args):
         retry_failed_runs=args['--retry_failed_runs'],
         ff=args['--ff'],
         implicit_water_model=args['--water_model'],
-        verbose=args['--verbose'],
+        loglevel=loglevel,
         **api_params
     )

ensembler/core.py

@@ -212,6 +212,11 @@ class ManualOverrides:
     """
     Reads in user-defined override data from a YAML file named "manual-overrides.yaml"
 
+    Parameters
+    ----------
+    manual_overrides_filepath: str
+        In normal use, this should not need to be set. Defaults to 'manual-overrides.yaml'
+
     Example file contents
     ---------------------
 
@@ -231,17 +236,27 @@ class ManualOverrides:
             - 4Q2A
             - 4CTB
             - 4QOX
+    refinement:
+        ph: 8.0
+        custom_residue_variants:
+            DDR1_HUMAN_D0_PROTONATED:
+                # keyed by 0-based residue index
+                35: ASH
+
+    Or see `ensembler/tests/example_project/manual_overrides.yaml` for an example file.
     """
-    def __init__(self):
-        import yaml
-        if os.path.exists(manual_overrides_filename):
-            with open(manual_overrides_filename, 'r') as manual_overrides_file:
+    def __init__(self, manual_overrides_filepath=None):
+        if not manual_overrides_filepath:
+            manual_overrides_filepath = manual_overrides_filename
+        if os.path.exists(manual_overrides_filepath):
+            with open(manual_overrides_filepath, 'r') as manual_overrides_file:
                 manual_overrides_yaml = yaml.load(manual_overrides_file, Loader=YamlLoader)
         else:
             manual_overrides_yaml = {}
 
         self.target = TargetManualOverrides(manual_overrides_yaml)
         self.template = TemplateManualOverrides(manual_overrides_yaml)
+        self.refinement = RefinementManualOverrides(manual_overrides_yaml)
 
 
 class TargetManualOverrides:
@@ -258,7 +273,7 @@ class TargetManualOverrides:
     """
     def __init__(self, manual_overrides_yaml):
         target_dict = manual_overrides_yaml.get('target-selection')
-        if target_dict != None:
+        if target_dict is not None:
             self.domain_spans = target_dict.get('domain-spans')
         else:
             self.domain_spans = {}
@@ -282,7 +297,7 @@ class TemplateManualOverrides:
     """
     def __init__(self, manual_overrides_yaml):
         template_dict = manual_overrides_yaml.get('template-selection')
-        if template_dict != None:
+        if template_dict is not None:
             self.min_domain_len = template_dict.get('min-domain-len')
             self.max_domain_len = template_dict.get('max-domain-len')
             self.domain_spans = template_dict.get('domain-spans')
@@ -294,6 +309,29 @@ def __init__(self, manual_overrides_yaml):
             self.skip_pdbs = []
 
 
+class RefinementManualOverrides:
+    """
+    Parameters
+    ----------
+    manual_overrides_yaml: dict
+
+    Attributes
+    ----------
+    ph: float or NoneType
+    custom_residue_variants_by_targetid: dict or NoneType
+        dict with structure {`targetid`: {residue_index: residue_name}, ...} where
+        e.g. {'DDR1_HUMAN_D0_PROTONATED': {35: 'ASH'}}
+    """
+    def __init__(self, manual_overrides_yaml):
+        refinement_dict = manual_overrides_yaml.get('refinement')
+        if refinement_dict is not None:
+            self.ph = refinement_dict.get('ph')
+            self.custom_residue_variants_by_targetid = refinement_dict.get('custom_residue_variants')
+        else:
+            self.ph = None
+            self.custom_residue_variants_by_targetid = {}
+
+
 def gen_metadata_filename(ensembler_stage, metadata_file_index):
     for modelling_stage in ['build_models', 'cluster_models', 'refine_implicit_md', 'solvate_models', 'determine_nwaters', 'refine_explicit_md']:
         if ensembler_stage == modelling_stage:

ensembler/initproject.py

@@ -3,20 +3,19 @@
 import sys
 import os
 import re
-
+import shutil
 from lxml import etree
 import Bio.SeqUtils
 import Bio.SeqIO
 from Bio.Seq import Seq
 from Bio.SeqRecord import SeqRecord
-
 import ensembler
 import ensembler.version
 import ensembler.targetexplorer
 import ensembler.uniprot
 import ensembler.pdb
-from ensembler.utils import file_exists_and_not_empty
-from ensembler.core import mpistate, logger
+from ensembler.utils import file_exists_and_not_empty, get_installed_resource_filename
+from ensembler.core import mpistate, logger, manual_overrides_filename
 
 
 class TemplateData:
@@ -43,6 +42,7 @@ def __init__(self, project_toplevel_dir, run_main=True):
     @ensembler.utils.notify_when_done
     def _init_project(self):
         self._create_project_dirs()
+        self._write_manual_overrides_file()
         self._write_init_metadata()
 
     def _create_project_dirs(self):
@@ -56,6 +56,16 @@ def _create_project_dirs(self):
         ensembler.utils.create_dir(os.path.join(self.project_toplevel_dir, ensembler.core.default_project_dirnames.templates_structures_resolved))
         ensembler.utils.create_dir(os.path.join(self.project_toplevel_dir, ensembler.core.default_project_dirnames.templates_structures_modeled_loops))
 
+    def _write_manual_overrides_file(self):
+        if not os.path.exists(manual_overrides_filename):
+            template_manual_overrides_filepath = get_installed_resource_filename(
+                os.path.join('resources', 'template-manual-overrides.yaml')
+            )
+            manual_overrides_filepath = os.path.join(
+                self.project_toplevel_dir, manual_overrides_filename
+            )
+            shutil.copy(template_manual_overrides_filepath, manual_overrides_filepath)
+
     def _write_init_metadata(self):
         project_metadata = ensembler.core.ProjectMetadata(project_stage='init', project_toplevel_dir=self.project_toplevel_dir)
         init_metadata = self._gen_init_metadata()

ensembler/modeling.py

@@ -135,7 +135,7 @@ def pdbfix_template(template_full_seq, overwrite_structures=False):
             template_full_seq.id + '.pdb'
         )
         fixer = pdbfixer.PDBFixer(filename=template_filepath)
-        chainid = fixer.structureChains[0].chain_id
+        chainid = next(fixer.structure.iter_chains()).chain_id
         seq_obj = simtk.openmm.app.internal.pdbstructure.Sequence(chainid)
         for r in template_full_seq.seq:
             resi3 = Bio.SeqUtils.seq3(r).upper()
@@ -179,6 +179,8 @@ def pdbfix_template(template_full_seq, overwrite_structures=False):
             'MPI rank %d pdbfixer error for template %s - see logfile' %
             (mpistate.rank, template_full_seq.id)
         )
+        logger.debug(e)
+        logger.debug(trbk)
 
 
 def remove_missing_residues_at_termini(fixer, len_full_seq):

ensembler/packaging.py

@@ -4,7 +4,7 @@
 from ensembler.core import get_targets_and_templates, select_templates_by_seqid_cutoff
 from ensembler.utils import set_loglevel, read_file_contents_gz_or_not
 import simtk.unit as unit
-import simtk.openmm as openmm
+import simtk.openmm as mm
 import mdtraj
 
 fah_projects_dir = os.path.join(default_project_dirnames.packaged_models, 'fah-projects')
@@ -15,9 +15,9 @@ def package_for_fah(process_only_these_targets=None,
                     template_seqid_cutoff=None,
                     nclones=1, archive=False,
                     openmm_platform='Reference',
-                    timestep=2.0 * unit.femtoseconds,
-                    collision_rate=1.0 / unit.picosecond,
                     temperature=300.0 * unit.kelvin,
+                    collision_rate=1.0 / unit.picosecond,
+                    timestep=2.0 * unit.femtoseconds,
                     loglevel=None):
     """
     Create the input files and directory structure necessary to start a Folding@Home project.
@@ -76,9 +76,9 @@ def package_for_fah(process_only_these_targets=None,
             system = setup_system_and_integrator_files(
                 target,
                 sorted_valid_templates[0],
-                timestep,
+                temperature,
                 collision_rate,
-                temperature
+                timestep
             )
 
             renumbered_resnums = get_renumbered_topol_resnums(target)
@@ -110,7 +110,7 @@ def package_for_fah(process_only_these_targets=None,
                 nclones,
                 temperature,
                 collision_rate,
-                temperature,
+                timestep,
                 openmm_platform,
                 renumbered_resnums,
             )
@@ -206,10 +206,11 @@ def create_target_project_dir(target):
 
 def setup_system_and_integrator_files(target,
                                       template,
-                                      timestep,
+                                      temperature,
                                       collision_rate,
-                                      temperature
+                                      timestep
                                       ):
+    logger.debug('Copying system and integrator files for template {}'.format(template.id))
     models_target_dir = os.path.join(default_project_dirnames.models, target.id)
     template_dir = os.path.join(models_target_dir, template.id)
     target_project_dir = os.path.join(fah_projects_dir, target.id)
@@ -218,10 +219,10 @@ def setup_system_and_integrator_files(target,
     dest_system_filepath = os.path.join(target_project_dir, 'system.xml')
     dest_integrator_filepath = os.path.join(target_project_dir, 'integrator.xml')
 
-    system = openmm.XmlSerializer.deserialize(
+    system = mm.XmlSerializer.deserialize(
         read_file_contents_gz_or_not(source_system_filepath)
     )
-    state = openmm.XmlSerializer.deserialize(
+    state = mm.XmlSerializer.deserialize(
         read_file_contents_gz_or_not(source_state_filepath)
     )
 
@@ -230,17 +231,17 @@ def setup_system_and_integrator_files(target,
     system.setDefaultPeriodicBoxVectors(*box_vectors)
 
     # Create new integrator to use.
-    integrator = openmm.LangevinIntegrator(temperature, collision_rate, timestep)
+    integrator = mm.LangevinIntegrator(temperature, collision_rate, timestep)
 
     # TODO: Make sure MonteCarloBarostat temperature matches set temperature.
 
     # Serialize System.
     with open(dest_system_filepath, 'w') as dest_system_file:
-        dest_system_file.write(openmm.XmlSerializer.serialize(system))
+        dest_system_file.write(mm.XmlSerializer.serialize(system))
 
     # Serialize Integrator
     with open(dest_integrator_filepath, 'w') as dest_integrator_file:
-        dest_integrator_file.write(openmm.XmlSerializer.serialize(integrator))
+        dest_integrator_file.write(mm.XmlSerializer.serialize(integrator))
 
     return system
 
@@ -337,7 +338,7 @@ def generate_fah_run(target_project_dir,
                     read_file_contents_gz_or_not(source_system_structure_filepath)
                 )
 
-        state = openmm.XmlSerializer.deserialize(
+        state = mm.XmlSerializer.deserialize(
             read_file_contents_gz_or_not(source_openmm_state_filepath)
         )
 
@@ -346,13 +347,12 @@ def generate_fah_run(target_project_dir,
             run_seqid_file.write(read_file_contents_gz_or_not(source_seqid_filepath))
 
         # Create new integrator to use.
-        integrator = openmm.LangevinIntegrator(temperature, collision_rate, timestep)
+        integrator = mm.LangevinIntegrator(temperature, collision_rate, timestep)
 
         # Create Context so we can randomize velocities.
-        platform = openmm.Platform.getPlatformByName(openmm_platform)
-        context = openmm.Context(system, integrator, platform)
+        platform = mm.Platform.getPlatformByName(openmm_platform)
+        context = mm.Context(system, integrator, platform)
         context.setPositions(state.getPositions())
-        context.setVelocities(state.getVelocities())
         box_vectors = state.getPeriodicBoxVectors()
         context.setPeriodicBoxVectors(*box_vectors)
 
@@ -371,7 +371,7 @@ def generate_fah_run(target_project_dir,
                 enforcePeriodicBox=True
             )
             with open(state_filename, 'w') as state_file:
-                state_file.write(openmm.XmlSerializer.serialize(state))
+                state_file.write(mm.XmlSerializer.serialize(state))
 
     except Exception as e:
         import traceback