Update MCSS calculation (#33)

* calculate RMSD using atom map instead of sub-Mols

* small bugfix

* bugfix when whole molecule is MCSS

Some issues occur with the RMSD calculation when the whole molecule is
the MCSS and map is set. The RMSD is different from what obrms reports.
Therefore, we keep track of whether or not the molecules are the same
and calculate the RMSD for these molecules without setting map.

* parallel mcss bugfixes

* assertion for atom_map length

* ridiculous parameter bug fix

* MCSS param available

* Benchmark Open-ComBind with new MCSS parameters

* docstring updates for MCSS updates

* Added Open-ComBind statistics

open_combind/cli/cli.py

@@ -12,7 +12,7 @@
 ###############################################################################
 
 # Defaults
-SHAPE_VERSION = 'pharm_max'
+MCSS_PARAM = 'strict'
 IFP_VERSION = 'rd1'
 
 
@@ -163,7 +163,7 @@ def dock_ligands(template, root, ligands, screen, slurm, now, dock_file):
 @click.argument('poseviewers', nargs=-1)
 @click.option('--native', default='structures/ligands/*_lig.sdf')
 @click.option('--ifp-version', default=IFP_VERSION)
-@click.option('--shape-version', default=SHAPE_VERSION)
+@click.option('--mcss-param', default=MCSS_PARAM)
 @click.option('--screen', is_flag=True)
 @click.option('--max-poses', default=100)
 @click.option('--no-mcss', is_flag=True)
@@ -173,12 +173,12 @@ def dock_ligands(template, root, ligands, screen, slurm, now, dock_file):
 @click.option('--template', default='structures/template/*.template')
 @click.option('--check-center-ligs', is_flag=True)
 def featurize(root, poseviewers, native, ifp_version,
-            shape_version, screen, no_mcss,
+            mcss_param, screen, no_mcss,
             use_shape, processes, max_poses, no_cnn, template, check_center_ligs):
 
     oc.featurize(root, poseviewers, native=native, no_mcss=no_mcss, use_shape=use_shape,
                 max_poses=max_poses, no_cnn=no_cnn, screen=screen, ifp_version=ifp_version,
-                shape_version=shape_version, processes=processes, template=template, check_center_ligs=check_center_ligs)
+                mcss_param=mcss_param, processes=processes, template=template, check_center_ligs=check_center_ligs)
 ################################################################################
 
 

open_combind/features/features.py

@@ -32,8 +32,8 @@ class Features:
         Path to root directory of poses/features
     ifp_version : str, default='rd1'
         Version of the interaction fingerprint to use for featurization
-    shape_version : str, default='pharm_max'
-        Version of the shape algorithm to use for featurization (not implemented)
+    mcss_param : str, default='strict'
+        MCSS parameter set to use for MCSS RMSD featurization, either ``'strict'`` or ``'relaxed'``
     max_poses : str, default=10000  
         Maximum number of poses per ligand used to compute features
     pv_root
@@ -47,8 +47,8 @@ class Features:
     ----------
     root : str
         Root directory where the features/poses should be looked for
-    ifp_version : str
-        Which version of the interaction fingerprint to use
+    mcss_param : str
+        Which parameters to use for MCSS featurization either ``'strict'`` or ``'relaxed'``
     shape_version : str
         (not currently used) which version of the shape algorithm to use for featurization
     max_poses : int
@@ -60,7 +60,7 @@ class Features:
     raw : dict
         Contains the raw data of the computed (or loaded) features where the key is the feature name
     """
-    def __init__(self, root, ifp_version='rd1', shape_version='pharm_max',
+    def __init__(self, root, ifp_version='rd1', mcss_param='strict',
                  max_poses=10000, pv_root=None,
                  ifp_features=['hbond', 'saltbridge', 'contact'], cnn_scores=True,
                  template='',check_center_ligs=False):
@@ -71,7 +71,7 @@ def __init__(self, root, ifp_version='rd1', shape_version='pharm_max',
             self.pv_root = pv_root
 
         self.ifp_version = ifp_version
-        self.shape_version = shape_version
+        self.mcss_param = mcss_param
         self.max_poses = max_poses
         self.ifp_features = ifp_features
         self.cnn_scores = cnn_scores
@@ -104,7 +104,7 @@ def get_molecules_from_files(self, pvs, native=False, center_ligand=None):
 
         Returns
         -------
-        molbundle_dict : :class:`dict[str, list[~rdkit.Chem.rdchem.Mol]]<dict>`
+        molbundle_dict : :class:`dict[str, list[Mol]]<dict>`
             Dictionary of pose viewer files to list of molecules
         """
         molbundle_dict = dict()
@@ -231,7 +231,7 @@ def load_single_features(self, pvs, ligands=None, center_ligand=None):
         ----------
         pvs : :class:`list[str]<list>`
             Poses that need features loaded
-        ligands : 
+        ligands : :class:`list[str]<list>`, default=None
 
 		center_ligand: :class:`~rdkit.Chem.rdchem.Mol`, default=None
 			Ligand to center the poses around
@@ -597,12 +597,12 @@ def compute_shape(self, poses1, poses2, out, processes=1):
         """
         if processes != 1:
             from open_combind.features.shape import shape_mp
-            sims = shape_mp(poses2, poses1, version=self.shape_version,processes=processes).T
+            sims = shape_mp(poses2, poses1,processes=processes).T
         else:
             from open_combind.features.shape import shape
             # More efficient to have longer pose list provided as second argument.
             # This only matters for screening.
-            sims = shape(poses2, poses1, version=self.shape_version).T
+            sims = shape(poses2, poses1).T
         np.save(out, sims)
 
     def compute_mcss(self, poses1, poses2, out, processes=1):
@@ -622,8 +622,8 @@ def compute_mcss(self, poses1, poses2, out, processes=1):
         """
         if processes != 1:
             from open_combind.features.mcss import mcss_mp
-            rmsds = mcss_mp(poses1, poses2, processes)
+            rmsds = mcss_mp(poses1, poses2, param_string=self.mcss_param, processes=processes)
         else:
             from open_combind.features.mcss import mcss
-            rmsds = mcss(poses1, poses2)
+            rmsds = mcss(poses1, poses2, param_string=self.mcss_param)
         np.save(out, rmsds)