ChemSpace

ChemSpace(data=None, data_indices=None, dask_cluster=None, dask_local_cluster_kwargs={})

Streamline working with many RMols or a specific chemical space by employing a pandas dataframe, in combination with Dask for parallellisation.

Methods:

• optimise_in_receptor –

  Return lists of energies.

• discard_missing –

  Remove from this list the molecules that have no conformers

• add_rgroups –

  Note that if they are Smiles:

• add_data –

  :param data: dictionary {"Smiles": [], "h": [], ... }

• add_smiles –

  Add a list of Smiles into this ChemicalSpace

• evaluate –

  :param indices:

• add_enamine_molecules –

  For the best scoring molecules, find similar molecules in Enamine REAL database

• active_learning –

  Model the data using the Training subset. Then use the active learning query method.

• compute_fps –

  :param smiles_tuple: It has to be a tuple to be hashable (to work with caching).

• to_sdf –

  Write every molecule and all its fields as properties, to an SDF file.

Source code in fegrow/package.py

def __init__(
    self,
    data=None,
    data_indices=None,
    dask_cluster=None,
    dask_local_cluster_kwargs={},
):
    if data is None:
        data = ChemSpace.DATAFRAME_DEFAULT_VALUES

    self.df = pandas.DataFrame(data, index=data_indices)

    ChemSpace._dask_cluster = dask_cluster

    if ChemSpace._dask_cluster is None:
        logger.info(
            "No Dask cluster configured. Creating a local cluster of threads. "
        )
        warnings.warn(
            "ANI uses TORCHAni which is not threadsafe, leading to random SEGFAULTS. "
            "Use a Dask cluster with processes as a work around "
            "(see the documentation for an example of this workaround) ."
        )

        kwargs = {
            "n_workers": None,
            "processes": False,  # turn off Nanny to avoid the problem
            # with loading of the main file (ie executing it)
            "dashboard_address": ":8989",
            **dask_local_cluster_kwargs,
        }
        ChemSpace._dask_cluster = LocalCluster(**kwargs)
        # ChemSpace._dask_cluster = Scheduler()
        # ChemSpace._dask_cluster = LocalCluster(preload_nanny=["print('Hi Nanny')"],
        #                                        preload=["pint"], n_workers=1
        #                                        ) #asynchronous=True)

    ChemSpace._dask_client = Client(
        ChemSpace._dask_cluster
    )  # ChemSpace._dask_cluster, asynchronous=True)
    print(f"Dask can be watched on {ChemSpace._dask_client.dashboard_link}")

    self._scaffolds = []
    self._model = None
    self._query = None
    self._query_label = None

optimise_in_receptor

optimise_in_receptor(*args, **kwargs)

Return lists of energies.

Source code in fegrow/package.py

def optimise_in_receptor(self, *args, **kwargs):
    """
    Return lists of energies.
    """

    # daskify parameters
    args = [dask.delayed(arg) for arg in args]
    kwargs = {k: dask.delayed(v) for k, v in kwargs.items()}

    # create the dask jobs
    delayed_optimise_in_receptor = dask.delayed(optimise_in_receptor)
    jobs = {}
    for i, row in self.df.iterrows():
        if row.Mol.GetNumConformers() == 0:
            print(
                f"Warning: mol {i} has no conformers. Ignoring receptor optimisation."
            )
            continue

        jobs[row.Mol] = delayed_optimise_in_receptor(row.Mol, *args, **kwargs)

    # dask batch compute
    results = dict(zip(jobs.keys(), self.dask_client.compute(list(jobs.values()))))

    # extract results
    dfs = []
    for mol, result in results.items():
        opt_mol, energies = result.result()
        mol.RemoveAllConformers()
        # replace the conformers with the optimised ones
        [mol.AddConformer(c) for c in opt_mol.GetConformers()]

        mol.SetProp("energies", str(energies))
        dfs.append(pandas.DataFrame({}))
        mol._save_opt_energies(energies)

discard_missing

discard_missing()

Remove from this list the molecules that have no conformers

Source code in fegrow/package.py

def discard_missing(self):
    """
    Remove from this list the molecules that have no conformers
    """
    ids_to_remove = []
    for i, row in self.df.iterrows():
        if row.Mol.GetNumConformers() == 0:
            print(f"Discarding a molecule (id {i}) due to the lack of conformers. ")
            ids_to_remove.append(i)

    self.df = self.df[~self.df.index.isin(ids_to_remove)]
    return ids_to_remove

add_rgroups

add_rgroups(rgroups_linkers, rgroups2=None, alltoall=False)

Note that if they are Smiles

• if they have an * atom (e.g. RDKit atom.SetAtomicNum(0)), this will be used for attachment to the scaffold
• if they don't have an * atom, the scaffold will be fitted as a substructure

First link the linker to the scaffold. Then add the rgroups.

:param rgroups2: A list of Smiles. Molecules will be accepted and converted to Smiles. :param linker: A molecule. Ideally it has 2 atatchement points. :return:

Source code in fegrow/package.py

def add_rgroups(self, rgroups_linkers, rgroups2=None, alltoall=False):
    """
    Note that if they are Smiles:
     - if they have an * atom (e.g. RDKit atom.SetAtomicNum(0)), this will be used for attachment to the scaffold
     - if they don't have an * atom, the scaffold will be fitted as a substructure

    First link the linker to the scaffold. Then add the rgroups.

    :param rgroups2: A list of Smiles. Molecules will be accepted and converted to Smiles.
    :param linker: A molecule. Ideally it has 2 atatchement points.
    :return:
    """
    scaffold = dask.delayed(self._scaffolds[0])

    if not isinstance(rgroups_linkers, typing.Iterable):
        rgroups_linkers = [rgroups_linkers]

    if rgroups2 is not None and not isinstance(rgroups2, typing.Iterable):
        rgroups2 = [rgroups2]

    # create the dask jobs
    delayed_build_molecule = dask.delayed(build_molecule)

    jobs = [delayed_build_molecule(scaffold, linker) for linker in rgroups_linkers]

    # if linkers and rgroups are attached, add them in two iterations
    if rgroups2 is not None and not alltoall:
        # for each attached linker, attach an rgroup with the same position
        jobs = [
            delayed_build_molecule(scaffold_linked, rgroup)
            for rgroup, scaffold_linked in itertools.zip_longest(
                rgroups2, jobs, fillvalue=jobs[0]
            )
        ]
    elif rgroups2 is not None and alltoall:
        jobs = [
            delayed_build_molecule(scaffold_linked, rgroup)
            for rgroup, scaffold_linked in itertools.product(rgroups2, jobs)
        ]

    results = self.dask_client.compute(jobs)
    built_mols = [r.result() for r in results]

    # get Smiles
    built_mols_smiles = [Chem.MolToSmiles(mol) for mol in built_mols]

    # extract the H indices used for attaching the scaffold
    hs = [mol.GetIntProp("attachment_point") for mol in built_mols]

    self.add_data({"Smiles": built_mols_smiles, "Mol": built_mols, "h": hs})

add_data

add_data(data)

:param data: dictionary {"Smiles": [], "h": [], ... } :return:

Source code in fegrow/package.py

def add_data(self, data):
    """

    :param data: dictionary {"Smiles": [], "h": [], ... }
    :return:
    """

    # ensure that the new indices start at the end
    last_index = max([self.df.index.max() + 1])
    if np.isnan(last_index):
        last_index = 0

    # ensure correct default values in the new rows
    data_with_defaults = ChemSpace.DATAFRAME_DEFAULT_VALUES.copy()
    data_with_defaults.update(data)

    # update the internal dataframe
    new_indices = range(last_index, last_index + len(data_with_defaults["Smiles"]))
    prepared_data = pandas.DataFrame(data_with_defaults, index=new_indices)
    self.df = pandas.concat([self.df, prepared_data])
    return prepared_data

add_smiles

add_smiles(smiles_list, h=NA, protonate=False)

Add a list of Smiles into this ChemicalSpace

:param h: which h was used to connect to the :param protonate: use openbabel to protonate each smile :return:

Source code in fegrow/package.py

def add_smiles(self, smiles_list, h=pandas.NA, protonate=False):
    """
    Add a list of Smiles into this ChemicalSpace

    :param h: which h was used to connect to the
    :param protonate: use openbabel to protonate each smile
    :return:
    """

    if protonate:
        delayed_protonations = [
            DaskTasks.obabel_protonate(smi) for smi in smiles_list
        ]
        jobs = self.dask_client.compute(delayed_protonations)
        smiles_list = [job.result() for job in jobs]

    # convert the Smiles into molecules
    params = Chem.SmilesParserParams()
    params.removeHs = False
    mols = [Chem.MolFromSmiles(smiles, params=params) for smiles in smiles_list]

    self.add_data({"Smiles": smiles_list, "Mol": mols, "h": h})

_evaluate_experimental

_evaluate_experimental(indices=None, num_conf=10, minimum_conf_rms=0.5, min_dst_allowed=1)

Generate the conformers and score the subset of molecules.

E.g. :param indices: The indices in the dataframe to be run through the pipeline. If None, all molecules are evaluated. :return:

Source code in fegrow/package.py

def _evaluate_experimental(
    self, indices=None, num_conf=10, minimum_conf_rms=0.5, min_dst_allowed=1
):
    """
    Generate the conformers and score the subset of molecules.

    E.g.
    :param indices: The indices in the dataframe to be run through the pipeline.
        If None, all molecules are evaluated.
    :return:
    """

    if len(self._scaffolds) == 0:
        print("Please add scaffolds to the system for the evaluation. ")
    elif len(self._scaffolds) > 1:
        raise NotImplementedError(
            "For now we only allow working with one scaffold. "
        )

    # carry out the full pipeline, generate conformers, etc.
    # note that we have to find a way to pass all the molecules
    # this means creating a pipeline of tasks for Dask,
    # a pipeline that is well conditional

    ## GENERATE CONFORMERS
    num_conf = dask.delayed(num_conf)
    minimum_conf_rms = dask.delayed(minimum_conf_rms)
    protein = dask.delayed(prody_package.parsePDB(self._protein_filename))
    protein_file = dask.delayed(self._protein_filename)
    min_dst_allowed = dask.delayed(min_dst_allowed)
    RMol._check_download_gnina()
    gnina_path = dask.delayed(os.path.join(RMol.gnina_dir, "gnina"))

    # functions
    delayed_generate_conformers = dask.delayed(generate_conformers)
    delayed_remove_clashing_confs = dask.delayed(RMol.remove_clashing_confs)
    delayed_gnina = dask.delayed(gnina)

    # create dask jobs
    jobs = {}
    for i, row in self.df.iterrows():
        generated_confs = delayed_generate_conformers(
            row.Mol, num_conf, minimum_conf_rms
        )
        removed_clashes = delayed_remove_clashing_confs(
            generated_confs, protein, min_dst_allowed=min_dst_allowed
        )
        jobs[i] = delayed_gnina(removed_clashes, protein_file, gnina_path)

    # run all jobs
    results = dict(zip(jobs.keys(), self.dask_client.compute(list(jobs.values()))))

    # gather the results
    for i, result in results.items():
        mol, cnnaffinities = result.result()

        # extract the conformers
        input_mol = self.df.Mol[i]
        input_mol.RemoveAllConformers()
        [input_mol.AddConformer(c) for c in mol.GetConformers()]

        # save the affinities so that one can retrace which conformer has the best energy
        input_mol.SetProp("cnnaffinities", str(cnnaffinities))

        self.df.score[i] = max(cnnaffinities)

    logger.info(f"Evaluated {len(results)} cases")

evaluate

evaluate(indices: Union[Sequence[int], DataFrame] = None, scoring_function=None, gnina_path=None, gnina_gpu=False, num_conf=50, minimum_conf_rms=0.5, penalty=NA, al_ignore_penalty=True, **kwargs)

:param indices: :param scoring_function: :param gnina_path: :param gnina_gpu: :param num_conf: :param minimum_conf_rms: :param penalty: :param al_ignore_penalty: :param kwargs: :return:

Source code in fegrow/package.py

def evaluate(
    self,
    indices: Union[Sequence[int], pandas.DataFrame] = None,
    scoring_function=None,
    gnina_path=None,
    gnina_gpu=False,
    num_conf=50,
    minimum_conf_rms=0.5,
    penalty=pd.NA,
    al_ignore_penalty=True,
    **kwargs,
):
    """

    :param indices:
    :param scoring_function:
    :param gnina_path:
    :param gnina_gpu:
    :param num_conf:
    :param minimum_conf_rms:
    :param penalty:
    :param al_ignore_penalty:
    :param kwargs:
    :return:
    """

    # evaluate all molecules if no indices are picked
    if indices is None:
        indices = slice(None)

    if isinstance(indices, pandas.DataFrame):
        if len(indices) <= 2:
            raise ValueError("Please provide at least 3 items")
        indices = indices.index

    selected_rows = self.df.loc[indices]

    # discard computed rows
    selected_rows = selected_rows[selected_rows.score.isna()]

    if len(self._scaffolds) == 0:
        print("Please add scaffolds to the system for the evaluation. ")
    elif len(self._scaffolds) > 1:
        raise NotImplementedError(
            "For now we only allow working with one scaffold. "
        )

    # should be enough to do it once, shared
    ## GENERATE CONFORMERS

    if gnina_path is not None:
        # gnina_path = delayed(os.path.join(RMol.gnina_dir, 'gnina'))
        RMol.set_gnina(os.path.join(RMol.gnina_dir, "gnina"))
    RMol._check_download_gnina()

    num_conf = dask.delayed(num_conf)
    minimum_conf_rms = dask.delayed(minimum_conf_rms)
    protein_file = dask.delayed(self._protein_filename)
    RMol._check_download_gnina()

    scaffold = dask.delayed(self._scaffolds[0])
    # extract which hydrogen was used for the attachement
    h_attachements = [
        a.GetIdx() for a in self._scaffolds[0].GetAtoms() if a.GetAtomicNum() == 0
    ]

    h_attachement_index = None
    if len(h_attachements) > 0:
        h_attachement_index = h_attachements[0]

    # create dask jobs
    delayed_evaluate = dask.delayed(_evaluate_atomic)
    jobs = {}
    for i, row in selected_rows.iterrows():
        jobs[i] = delayed_evaluate(
            scaffold,
            row.Smiles,
            protein_file,
            h=h_attachement_index,
            num_conf=num_conf,
            minimum_conf_rms=minimum_conf_rms,
            scoring_function=scoring_function,
            gnina_gpu=gnina_gpu,
            **kwargs,
        )

    # run all
    results = dict(zip(jobs.keys(), self.dask_client.compute(list(jobs.values()))))

    # gather the results
    for i, result in results.items():
        Training = True
        build_succeeded = True

        try:
            mol, data = result.result()

            # save all data generated
            if mol is not None:
                for k, v in data.items():
                    mol.SetProp(k, str(v))

                # replace the original molecule with the new one
                self.df.at[i, "Mol"] = mol

            # extract the score
            score = data["score"]
        except subprocess.CalledProcessError as E:
            logger.error("Failed Process", E, E.cmd, E.output, E.stdout, E.stderr)
            score = penalty
            build_succeeded = False

            if al_ignore_penalty:
                Training = False
        except Exception:
            # failed to finish the protocol, set the penalty
            score = penalty
            build_succeeded = False

            if al_ignore_penalty:
                Training = False

        self.df.loc[i, ["score", "Training", "Success"]] = (
            score,
            Training,
            build_succeeded,
        )

    logger.info(f"Evaluated {len(results)} cases")
    return self.df.loc[indices]

add_enamine_molecules

add_enamine_molecules(n_best=1, results_per_search=100, remove_scaffold_h=False)

For the best scoring molecules, find similar molecules in Enamine REAL database and add them to the dataset.

Make sure you have the permission/license to use https://sw.docking.org/search.html this way.

@scaffold: The scaffold molecule that has to be present in the found molecules. If None, this requirement will be ignored. @molecules_per_smile: How many top results (molecules) per Smiles searched.

Source code in fegrow/package.py

def add_enamine_molecules(
    self, n_best=1, results_per_search=100, remove_scaffold_h=False
):
    """
    For the best scoring molecules, find similar molecules in Enamine REAL database
     and add them to the dataset.

    Make sure you have the permission/license to use https://sw.docking.org/search.html
        this way.

    @scaffold: The scaffold molecule that has to be present in the found molecules.
        If None, this requirement will be ignored.
    @molecules_per_smile: How many top results (molecules) per Smiles searched.
    """

    from pydockingorg import Enamine

    if len(self._scaffolds) > 1:
        raise NotImplementedError("Only one scaffold is supported atm.")

    scaffold = self._scaffolds[0]

    # get the best performing molecules
    vl = self.df.sort_values(by="score", ascending=False)
    best_vl_for_searching = vl[:n_best]

    # nothing to search for yet
    if len(best_vl_for_searching[~best_vl_for_searching.score.isna()]) == 0:
        return

    if len(set(best_vl_for_searching.h)) > 1:
        raise NotImplementedError("Multiple growth vectors are used. ")

    # filter out previously queried molecules
    new_searches = best_vl_for_searching[
        best_vl_for_searching.enamine_searched == False  # noqa: E712
    ]
    smiles_to_search = list(new_searches.Smiles)

    start = time.time()
    print(f"Querying Enamine REAL. Looking up {len(smiles_to_search)} smiles.")
    try:
        with Enamine() as DB:
            results: pandas.DataFrame = DB.search_smiles(
                smiles_to_search,
                remove_duplicates=True,
                results_per_search=results_per_search,
            )
    except requests.exceptions.HTTPError as HTTPError:
        print("Enamine API call failed. ", HTTPError)
        return
    print(
        f"Enamine returned with {len(results)} rows in {time.time() - start:.1f}s."
    )

    # update the database that this molecule has been searched
    self.df.loc[new_searches.index, "enamine_searched"] = True

    if len(results) == 0:
        print("The server did not return a single Smiles!")
        return

    # prepare the scaffold for testing its presence
    # specifically, the hydrogen was replaced and has to be removed
    # for now we assume we only are growing one vector at a time - fixme
    if remove_scaffold_h:
        scaffold_noh = Chem.EditableMol(scaffold)
        scaffold_noh.RemoveAtom(int(best_vl_for_searching.iloc[0].h))
        scaffold = scaffold_noh.GetMol()

    dask_scaffold = dask.delayed(scaffold)

    start = time.time()
    # protonate and check for scaffold
    delayed_protonations = [
        DaskTasks.obabel_protonate(smi.rsplit(maxsplit=1)[0])
        for smi in results.hitSmiles.values
    ]
    jobs = self.dask_client.compute(
        [
            DaskTasks.scaffold_check(smih, dask_scaffold)
            for smih in delayed_protonations
        ]
    )
    scaffold_test_results = [job.result() for job in jobs]
    scaffold_mask = [r[0] for r in scaffold_test_results]
    # smiles None means that the molecule did not have our scaffold
    protonated_smiles = [r[1] for r in scaffold_test_results if r[1] is not None]
    print(
        f"Dask obabel protonation + scaffold test finished in {time.time() - start:.2f}s."
    )
    print(
        f"Tested scaffold presence. Kept {sum(scaffold_mask)}/{len(scaffold_mask)}."
    )

    if len(scaffold_mask) > 0:
        similar = results[scaffold_mask]
        similar.hitSmiles = protonated_smiles
    else:
        similar = pandas.DataFrame(columns=results.columns)

    # filter out Enamine molecules which were previously added
    new_enamines = similar[~similar.id.isin(vl.enamine_id)]

    warnings.warn(
        f"Only one H vector is assumed and used. Picking {vl.h[0]} hydrogen on the scaffold. "
    )
    new_data = {
        "Smiles": list(new_enamines.hitSmiles.values),
        "h": vl.h[0],  # fixme: for now assume that only one vector is used
        "enamine_id": list(new_enamines.id.values),
    }

    print("Adding: ", len(new_enamines.hitSmiles.values))
    return self.add_data(new_data)

active_learning

active_learning(n=1, first_random=True, score_higher_better=True, model=None, query=None, learner_type=None)

Model the data using the Training subset. Then use the active learning query method.

See properties "model" and "query" for finer control.

It's better to save the FPs in the dataframe. Or in the underlying system. :return:

Source code in fegrow/package.py

def active_learning(
    self,
    n=1,
    first_random=True,
    score_higher_better=True,
    model=None,
    query=None,
    learner_type=None,
):
    """
    Model the data using the Training subset. Then use the active learning query method.

    See properties "model" and "query" for finer control.

    It's better to save the FPs in the dataframe. Or in the underlying system.
    :return:
    """

    training = self.df[self.df.Training]
    selection = self.df[~self.df.Training]

    if training.empty:
        if first_random:
            warnings.warn(
                "Selecting randomly the first samples to be studied (no score data yet). "
            )
            return selection.sample(n)
        else:
            raise ValueError(
                'There is no scores for active learning. Please use the "first_random" property. '
            )

    # get the scored subset
    # fixme - multitarget?
    train_targets = training["score"].to_numpy(dtype=float)

    library_features = self.compute_fps(tuple(self.df.Smiles))

    train_features = library_features[training.index]

    selection_features = library_features[selection.index]

    import fegrow.al

    if model is not None:
        self.model = model
    if self.model is None:
        self.model = fegrow.al.Model.gaussian_process()

    if query is not None:
        self.query = query

    # employ Greedy query by default
    if self.query is None:
        self.query = fegrow.al.Query.Greedy()

    # update on how many to querry
    query = functools.partial(self.query, n_instances=n)

    target_multiplier = 1
    if score_higher_better is True:
        target_multiplier = -1

    if self.query_label in ["greedy", "thompson", "EI", "PI"]:
        target_multiplier *= 1
    elif self.query_label == "UCB":
        target_multiplier *= -1

    train_targets = train_targets * target_multiplier

    # only GP uses Bayesian Optimizer
    if learner_type is not None:
        learner = learner_type(
            estimator=self.model,
            X_training=train_features,
            y_training=train_targets,
            query_strategy=query,
        )
    elif isinstance(self.model, gaussian_process.GaussianProcessRegressor):
        learner = modAL.models.BayesianOptimizer(
            estimator=self.model,
            X_training=train_features,
            y_training=train_targets,
            query_strategy=query,
        )
    else:
        learner = modAL.models.ActiveLearner(
            estimator=self.model,
            X_training=train_features,
            y_training=train_targets,
            query_strategy=query,
        )

    inference = learner.predict(library_features) * target_multiplier

    self.df["regression"] = inference.T.tolist()

    selection_idx, _ = learner.query(selection_features)

    return selection.iloc[selection_idx]

compute_fps cached

compute_fps(smiles_tuple)

:param smiles_tuple: It has to be a tuple to be hashable (to work with caching). :return:

Source code in fegrow/package.py

@functools.cache
def compute_fps(self, smiles_tuple):
    """
    :param smiles_tuple: It has to be a tuple to be hashable (to work with caching).
    :return:
    """
    futures = self._dask_client.map(ChemSpace._compute_fp_from_smiles, smiles_tuple)
    fps = np.array([r.result() for r in futures])

    return fps

to_sdf

to_sdf(filename, failed=False, unbuilt=True)

Write every molecule and all its fields as properties, to an SDF file.

:return:

Source code in fegrow/package.py

def to_sdf(self, filename, failed=False, unbuilt=True):
    """
    Write every molecule and all its fields as properties, to an SDF file.

    :return:
    """
    with Chem.SDWriter(filename) as SD:
        columns = self.df.columns.to_list()
        columns.remove("Mol")

        for i, row in self.df.iterrows():
            # ignore this molecule because it failed during the build
            if failed is False and row.Success is False:
                continue

            # ignore this molecule because it was not built yet
            if unbuilt is False and row.Success is pandas.NA:
                continue

            mol = row.Mol
            mol.SetIntProp("index", i)
            for column in columns:
                value = getattr(row, column)
                mol.SetProp(column, str(value))

            mol.ClearProp("attachement_point")
            SD.write(mol)