Merge pull request #76 from FilomenoSanchez/master

Adding AMPLE helical ensembles

ample/main.py

@@ -242,6 +242,12 @@ def ensembling(self, optd):
         elif optd['ideal_helices']:
             ample_util.ideal_helices(optd)
             logger.info("*** Using ideal helices to solve structure ***")
+            logger.warning('If ideal helices do not solve the structure, you may want to use -helical_ensembles in '
+                           'place of -ideal_helices. AMPLE will then use a new set of helical ensembles which has been '
+                           'very successful on solving challenging cases!')
+        elif optd['helical_ensembles']:
+            ample_util.ideal_helices(optd, ensembles=True)
+            logger.info("*** Using helical ensembles to solve structure ***")
         else:
             # Check we have some models to work with
             if not (optd['single_model_mode'] or optd['processed_models']):
@@ -415,8 +421,8 @@ def monitor():
 
         if not ok:
             msg = (
-                "An error code was returned after running MRBUMP on the ensembles!\n"
-                + "For further information check the logs in directory: {0}".format(optd['mrbump_dir'])
+                    "An error code was returned after running MRBUMP on the ensembles!\n"
+                    + "For further information check the logs in directory: {0}".format(optd['mrbump_dir'])
             )
             logger.critical(msg)
 

ample/util/ample_util.py

@@ -266,14 +266,17 @@ def filename_append(filename=None, astr=None, directory=None, separator="_"):
     return os.path.join(directory, name)
 
 
-def ideal_helices(optd):
+def ideal_helices(optd, ensembles=False):
     """Get some ideal helices
 
     Parameters
     ----------
     nresidues : int
        Number of residues to be used
 
+    ensembles : bool
+       If True then helical ensembles will be used instead of polyalanine ideal helices
+
     Returns
     -------
     pdbs : list
@@ -282,13 +285,50 @@ def ideal_helices(optd):
     """
     nresidues = optd['fasta_length']
     include_dir = os.path.join(SHARE_DIR, 'include')
-    names = ['polyala5', 'polyala10', 'polyala15', 'polyala20', 'polyala25', 'polyala30', 'polyala35', 'polyala40']
-    polya_lengths = [5, 10, 15, 20, 25, 30, 35, 40]
+    if not ensembles:
+        names = ['polyala_5', 'polyala_10', 'polyala_15', 'polyala_20', 'polyala_25', 'polyala_30', 'polyala_35',
+                 'polyala_40']
+    else:
+        names = ['ensemble_20_bfactor2_homogenous', 'ensemble_20_bfactor2_heterogenous',
+                 'ensemble_20_bfactor1_heterogenous', 'ensemble_20_bfactor1_homogenous',
+                 'ensemble_20_bfactor4_heterogenous', 'ensemble_20_bfactor3_homogenous',
+                 'ensemble_20_bfactor3_heterogenous', 'ensemble_20_bfactor4_homogenous',
+                 'ensemble_15_bfactor2_heterogenous', 'ensemble_15_bfactor2_homogenous',
+                 'ensemble_15_bfactor1_heterogenous', 'ensemble_15_bfactor1_homogenous',
+                 'ensemble_15_bfactor4_heterogenous',
+                 'ensemble_15_bfactor3_heterogenous', 'ensemble_15_bfactor4_homogenous',
+                 'ensemble_15_bfactor3_homogenous', 'ensemble_25_bfactor2_homogenous',
+                 'ensemble_25_bfactor2_heterogenous', 'ensemble_25_bfactor1_heterogenous',
+                 'ensemble_25_bfactor1_homogenous', 'ensemble_25_bfactor4_homogenous',
+                 'ensemble_25_bfactor3_homogenous', 'ensemble_25_bfactor4_heterogenous',
+                 'ensemble_25_bfactor3_heterogenous', 'ensemble_30_bfactor2_homogenous',
+                 'ensemble_30_bfactor2_heterogenous', 'ensemble_30_bfactor1_heterogenous',
+                 'ensemble_30_bfactor1_homogenous', 'ensemble_30_bfactor4_homogenous',
+                 'ensemble_30_bfactor4_heterogenous',
+                 'ensemble_30_bfactor3_heterogenous', 'ensemble_30_bfactor3_homogenous',
+                 'ensemble_35_bfactor2_homogenous', 'ensemble_35_bfactor2_heterogenous',
+                 'ensemble_35_bfactor1_heterogenous', 'ensemble_35_bfactor1_homogenous',
+                 'ensemble_35_bfactor4_heterogenous',
+                 'ensemble_35_bfactor4_homogenous', 'ensemble_35_bfactor3_homogenous',
+                 'ensemble_35_bfactor3_heterogenous', 'ensemble_40_bfactor2_homogenous',
+                 'ensemble_40_bfactor2_heterogenous', 'ensemble_40_bfactor1_heterogenous',
+                 'ensemble_40_bfactor1_homogenous', 'ensemble_40_bfactor4_heterogenous',
+                 'ensemble_40_bfactor3_heterogenous', 'ensemble_40_bfactor4_homogenous',
+                 'ensemble_40_bfactor3_homogenous', 'ensemble_10_bfactor2_heterogenous',
+                 'ensemble_10_bfactor2_homogenous', 'ensemble_10_bfactor1_heterogenous',
+                 'ensemble_10_bfactor1_homogenous', 'ensemble_10_bfactor3_homogenous',
+                 'ensemble_10_bfactor4_homogenous', 'ensemble_10_bfactor3_heterogenous',
+                 'ensemble_10_bfactor4_heterogenous', 'ensemble_5_bfactor2_heterogenous',
+                 'ensemble_5_bfactor2_homogenous', 'ensemble_5_bfactor1_heterogenous',
+                 'ensemble_5_bfactor1_homogenous', 'ensemble_5_bfactor4_homogenous',
+                 'ensemble_5_bfactor3_heterogenous', 'ensemble_5_bfactor3_homogenous',
+                 'ensemble_5_bfactor4_heterogenous']
 
     ensemble_options = {}
     ensembles_data = []
     pdbs = []
-    for name, nres in zip(names, polya_lengths):
+    for name in names:
+        nres = int(name.split('_')[1])
         ncopies = nresidues / nres
         if ncopies < 1:
             ncopies = 1

ample/util/argparse_util.py

@@ -152,6 +152,14 @@ def add_general_options(parser=None):
         metavar='True/False',
         help='Use ideal polyalanine helices to solve structure (8 helices: from 5-40 residues)',
     )
+    parser.add_argument(
+        '-helical_ensembles',
+        action=BoolAction,
+        nargs='?',
+        metavar='True/False',
+        help='Use helical ensembles to solve structure (64 ensembles)',
+    )
+
     parser.add_argument(
         '-improve_template', metavar='improve_template', help='Path to a template to improve - NMR, homolog'
     )

ample/util/benchmark_util.py

@@ -12,7 +12,8 @@
 import shutil
 import sys
 
-from ample.util import ample_util, csymmatch, mtz_util, pdb_edit, pdb_model, reforigin, residue_map, rio, shelxe, tm_util
+from ample.util import ample_util, csymmatch, mtz_util, pdb_edit, pdb_model, reforigin, residue_map, rio, shelxe, \
+    tm_util
 
 logger = logging.getLogger(__name__)
 
@@ -140,9 +141,8 @@ def analyse(amoptd, newroot=None):
         mtz_util.to_hkl(amoptd['mtz'], hkl_file=os.path.join(amoptd['benchmark_dir'], SHELXE_STEM + ".hkl"))
         shutil.copyfile(amoptd['native_pdb_std'], os.path.join(amoptd['benchmark_dir'], SHELXE_STEM + ".ent"))
 
-    if amoptd['native_pdb'] and not (
-        amoptd['homologs'] or amoptd['ideal_helices'] or amoptd['import_ensembles'] or amoptd['single_model_mode']
-    ):
+    if amoptd['native_pdb'] and not (amoptd['homologs'] or amoptd['ideal_helices'] or amoptd['helical_ensembles']
+                                     or amoptd['import_ensembles'] or amoptd['single_model_mode']):
         analyseModels(amoptd)
 
     # Get the ensembling data
@@ -176,7 +176,7 @@ def analyse(amoptd, newroot=None):
 
         # Hack for ideal helices where num_residues are missing
         if amoptd['ideal_helices'] and ('num_residues' not in d or d['num_residues'] is None):
-            d['num_residues'] = int(d['ensemble_name'].lstrip('polyala'))
+            d['num_residues'] = int(d['ensemble_name'].lstrip('polyala_'))
 
         # Get the ensemble data and add to the MRBUMP data
         d['ensemble_percent_model'] = int((float(d['num_residues']) / float(amoptd['fasta_length'])) * 100)
@@ -233,7 +233,6 @@ def analyse(amoptd, newroot=None):
 
 
 def analyseModels(amoptd):
-
     # Get hold of a full model so we can do the mapping of residues
     refModelPdb = glob.glob(os.path.join(amoptd['models_dir'], "*.pdb"))[0]
 
@@ -333,7 +332,6 @@ def analysePdb(amoptd):
 
 
 def analyseSolution(amoptd, d, mrinfo):
-
     logger.info("Benchmark: analysing result: {0}".format(d['ensemble_name']))
 
     mrPdb = None
@@ -395,11 +393,12 @@ def analyseSolution(amoptd, d, mrinfo):
         # to compare to the native to allow us to determine which parts of the ensemble correspond to which parts of
         # the native structure - or if we were unable to calculate a res_seq_map
         if not (
-            amoptd['homologs']
-            or amoptd['ideal_helices']
-            or amoptd['import_ensembles']
-            or amoptd['single_model_mode']
-            or amoptd['res_seq_map']
+                amoptd['homologs']
+                or amoptd['ideal_helices']
+                or amoptd['helical_ensembles']
+                or amoptd['import_ensembles']
+                or amoptd['single_model_mode']
+                or amoptd['res_seq_map']
         ):
 
             # Get reforigin info

ample/util/mrbump_util.py

@@ -555,7 +555,7 @@ def finalSummary(amoptd):
         return "Could not find any MRBUMP results in directory: {0}!".format(amoptd['mrbump_dir'])
 
     if 'ensembles_data' in amoptd and not (
-        amoptd['ideal_helices'] or amoptd['homologs'] or amoptd['single_model_mode']
+        amoptd['ideal_helices'] or amoptd['helical_ensembles'] or amoptd['homologs'] or amoptd['single_model_mode']
     ):
         results = []
         # Merge dictionaries together

ample/util/options_processor.py

@@ -258,7 +258,7 @@ def process_modelling_options(optd):
         optd['models'] = optd['cluster_dir']
         optd['make_frags'] = False
         optd['make_models'] = False
-    elif optd['ideal_helices']:
+    elif optd['ideal_helices'] or optd['helical_ensembles']:
         optd['make_frags'] = False
         optd['make_models'] = False
     elif optd['homologs']:

ample/util/pyrvapi_results.py

@@ -470,7 +470,8 @@ def _create_summary_tab(self):
 
     def do_create_ensembles_section(self, ample_dict):
         return (
-            not (ample_dict.get('single_model_mode') or ample_dict.get('homologs') or ample_dict.get('ideal_helices'))
+            not (ample_dict.get('single_model_mode') or ample_dict.get('homologs') or
+                 ample_dict.get('ideal_helices') or ample_dict.get('helical_ensembles'))
             and bool(ample_dict.get('ensembles_data'))
             and not self.summary_tab_ensemble_sec_id
         )

ample/util/reference_manager.py

@@ -79,7 +79,7 @@ def setup_sections(self, optd):
                 self.section_labels[section] = ['AMPLE', 'CCP4', 'AMPLE_COILED-COILS', 'AMPLE_CONTACTS']
             elif section == self.SECTIONS.MODELLING:
                 labels = []
-                if optd.get('ideal_helices'):
+                if optd.get('ideal_helices') or optd.get('helical_ensembles'):
                     labels.append('AMPLE_COILED-COILS')
                 if optd.get('make_models'):
                     labels.append('ROSETTA')

docs/examples/rst/ideal_helices.rst

@@ -50,6 +50,10 @@ Finally we provide some options about how AMPLE will run:
 
 For a full list possible options see :ref:`AMPLE options <cl_options>`.
 
+.. note::
+   If ideal helices do not solve the structure, you may want to use ``-helical_ensembles`` in place of ``-ideal_helices``. AMPLE will then use a new set of helical ensembles which has been very successful on solving challenging cases!
+
+
 ------------------------------------------------------------------
 
 AMPLE Output