Refactor kinmodel set of scripts

kinoml/modeling/complexes.py

@@ -0,0 +1,18 @@
+def create_protein_ligand_complex(protein, ligand):
+    """
+    Refactor this function to minimize IO
+
+    >>> def protein_ligand_concatenation (template_hits, target_fasta_file, ligands):
+    >>>     templates_df = pd.read_csv(template_hits, sep=",")
+    >>>     # first templatePDB in a template Column
+    >>>     top_first_hit_pdb = templates_df["template"].iloc[0]
+    >>>     target_pdb_path = os.path.dirname(template_hits)
+    >>>     target_pdb_name = os.path.basename(target_fasta_file).split('.')[0]
+    >>>     top_hit_model = '{0}_{1}.pdb'.format(target_pdb_name, top_first_hit_pdb)
+    >>>     for lig in ligands:
+    >>>         protein_model = os.path.join(target_pdb_path, top_hit_model)
+    >>>         basename_ligand = os.path.splitext(os.path.basename(lig))[0]
+    >>>         complex_protein_ligand = '{0}_{1}.pdb'.format(top_hit_model.split('.')[0], basename_ligand)
+    >>>         print('cat {0} {1} > protein_ligand_complex_top_1_comp_model/{2}'.format(protein_model, lig, complex_protein_ligand))
+    >>>         os.system('cat {0} {1} > protein_ligand_complex_top_1_comp_model/{2}'.format(protein_model, lig, complex_protein_ligand))
+    """

kinoml/modeling/ligands.py

@@ -0,0 +1,64 @@
+"""
+Modeling tools for small compounds
+"""
+
+import os
+from openforcefield.topology import Molecule
+from openforcefield.utils.toolkits import OpenEyeToolkitWrapper
+
+
+def load_molecule(file_or_smiles):
+    if os.path.isfile(file_or_smiles):
+        return Molecule.from_file(file_or_smiles)
+    return Molecule.from_smiles(file_or_smiles)
+
+
+def generate_conformers(molecule):
+    """
+    Refactor this so it uses the Python API
+
+    >>> def conf_gen(smi,maxconfs=2000):
+    >>>     smi_prefix = os.path.splitext(os.path.basename(smi))[0]
+    >>>     print('{0} -in {1}/{2} -out {1}/OMEGA/{3}_omega.sdf -prefix {1}/OMEGA/{3}_omega -warts true -maxconfs {4} -strict false'.format(OMEGA, os.getcwd(), smi, smi_prefix, maxconfs))
+    >>>     os.system('{0} -in {1}/{2} -out {1}/OMEGA/{3}_omega.sdf -prefix {1}/OMEGA/{3}_omega -warts true -maxconfs {4} -strict false'.format(OMEGA, os.getcwd(), smi, smi_prefix, maxconfs))
+    """
+    # read here: https://docs.eyesopen.com/toolkits/python/omegatk/omegaexamples.html
+    # or simply use the Openforcefield wrapper
+
+
+def superpose(molecule, *targets):
+    """
+    Refactor this so it uses the Python API
+
+    >>> def lig_alignment(conformer, template_database, rocs_maxconfs_output=100):
+    >>>     sdf_prefix = os.path.basename(os.path.splitext(conformer)[0]).split('_')[0]
+    >>>     for template in template_database:
+    >>>         template_id = "_".join(os.path.basename(template).split("_")[0:3])
+    >>>         print('{0} -dbase {1}/{2} -query {3} -prefix {4}_{5}_rocs -oformat sdf -maxconfs 30 -outputquery false -qconflabel title -outputdir {1}/ROCS/'.format(ROCS, os.getcwd(),conformer, template, sdf_prefix, template_id))
+    >>>         os.system('{0} -dbase {1}/{2} -query {3} -prefix {4}_{5}_rocs -oformat sdf -maxconfs 30 -outputquery false -qconflabel title -outputdir {1}/ROCS/'.format(ROCS, os.getcwd(),conformer, template, sdf_prefix, template_id))
+    """
+    if not OpenEyeToolkitWrapper.is_available():
+        raise RuntimeError("OpenEye Toolkit must be installed and licensed")
+    from openeye import oeshape
+    # read here: https://docs.eyesopen.com/toolkits/python/shapetk/shape_examples.html#rocs
+
+
+def parameterize_for_rosetta(molecule):
+    """
+    Refactor this so it uses the Python API
+
+    >>> def sdftoparams(mol2params, top_hits_sdf_path):
+    >>>     for file in top_hits_sdf_path:
+    >>>         out_put_file_name = os.path.splitext(os.path.basename(file))[0]
+    >>>         print('{0} {1} -p sdf2params/{2}'.format(mol2params, file, out_put_file_name))
+    >>>         os.system('{0} {1} -p sdf2params/{2}'.format(mol2params, file, out_put_file_name))
+    """
+    # import importlib.util
+    # from distutils.spawn import find_executable
+    # spec = importlib.util.spec_from_file_location("molfile_to_params", find_executable("molfile_to_params.py"))
+    # molfile_to_params = importlib.util.module_from_spec(spec)
+    # spec.loader.exec_module(molfile_to_params)
+
+    # from molfile_to_params import main as generate_rosetta_params
+    # args = []
+    # generate_rosetta_params(args)

kinoml/modeling/receptors.py

@@ -0,0 +1,93 @@
+"""
+Modeling tools for macromolecular receptors, usually proteins
+"""
+
+
+def align_sequence(sequence, *targets):
+    """
+    Refactor this (only available as an executable; in the future we might
+    want to have a `kinoml.utils.CommandlineWrapper` thingy)
+
+    >>> def emboss_needle_search(emboss_needle, target_seq_path, template_seq_path):
+    >>>     for template_seq in template_seq_path:
+    >>>         target_seq_id = os.path.basename(target_seq_path).split('.')[0]
+    >>>         template_seq_id = os.path.basename(template_seq).split('.')[0]
+    >>>         print('{0} -sid1 {1} -asequence {2}/{3} -sid2 {4} -bsequence {5} -gapopen 10.0 -gapextend 0.5 -aformat3 markx3 -outfile {2}/protein_comp_modeling/protein_seq_alignment_files/{1}_{4}.needle'.format(emboss_needle, target_seq_id, os.getcwd(), target_seq_path, template_seq_id, template_seq))
+    >>>         os.system('{0} -sid1 {1} -asequence {2}/{3} -sid2 {4} -bsequence {5} -gapopen 10.0 -gapextend 0.5 -aformat3 markx3 -outfile {2}/protein_comp_modeling/protein_seq_alignment_files/{1}_{4}.needle'.format(emboss_needle, target_seq_id, os.getcwd(), target_seq_path, template_seq_id, template_seq))
+
+    """
+
+
+def build_model_with_rosetta(query_sequence, alignment, structure_templates):
+    """
+    Refactor this to minimize file IO
+
+    >>> templates = pd.read_csv(template_hits, sep=",")
+    >>> target_seq = os.path.basename(target_seq_path).split('.')[0]
+    >>> templates['tar_tem_seq_alin'] = templates['template'].apply(lambda x: "{}_{}.needle".format(target_seq, x))
+    >>> templates['tar_tem_seq_alin'] = alignment_file_path+templates['tar_tem_seq_alin']
+    >>> top_hit_template_file_path = templates['tar_tem_seq_alin'].tolist()
+    >>>
+    >>> aligned_seq = defaultdict(list)
+    >>> for path in top_hit_template_file_path:
+    >>>     target_template_file_name = os.path.splitext(os.path.basename(path))[0]
+    >>>     target_name_fasta_format = '>{} ..'.format(target_template_file_name.split('_')[0])
+    >>>     template_name_fasta_format = '>{} ..'.format('_'.join(target_template_file_name.split('_')[1:]))
+    >>>     target_aligned_seq = ''
+    >>>     template_aligned_seq = ''
+    >>>     with open (path, 'r') as readFile:
+    >>>         parse = False
+    >>>         parse2 = False
+    >>>         for line in readFile:
+    >>>             line = line.strip()
+    >>>             if not parse:
+    >>>                 if line.startswith(target_name_fasta_format):
+    >>>                     parse = True
+    >>>             elif line.startswith(template_name_fasta_format):
+    >>>                 parse = False
+    >>>             else:
+    >>>                 target_aligned_seq+=line
+    >>>
+    >>>             if not parse2:
+    >>>                 if line.startswith(template_name_fasta_format):
+    >>>                     parse2 = True
+    >>>             elif line.startswith('#'):
+    >>>                 parse2 = False
+    >>>             else:
+    >>>                 template_aligned_seq += line
+    >>>     aligned_seq[target_template_file_name].append(target_aligned_seq)
+    >>>     aligned_seq[target_template_file_name].append(template_aligned_seq)
+    >>>
+    >>> target_seq_for_modeling = {}
+    >>> for name, alignment_file in aligned_seq.items():
+    >>>     top_hits_alignment = '{}\n{}\n{}\n\n'.format(name, alignment_file[0], alignment_file[1])
+    >>>     with open('protein_comp_modeling/top_hits_alignment.txt', 'a') as writeFile:
+    >>>         writeFile.write(top_hits_alignment)
+    >>>     target_seq_based_on_temp_pdb = ''
+    >>>     for i in range(len(alignment_file[0])):
+    >>>         if not alignment_file[1][i] == '-':
+    >>>             target_seq_based_on_temp_pdb += alignment_file[0][i]
+    >>>     target_seq_for_modeling[name]=target_seq_based_on_temp_pdb
+    >>>
+    >>> final_target_template_for_modeling = {}
+    >>> for target_template, target_final_seq in target_seq_for_modeling.items():
+    >>>     template_name = '_'.join(target_template.split('_')[1:])
+    >>>     temp_list_dir = os.listdir(template_pdb_path)
+    >>>     for template_hit in temp_list_dir:
+    >>>         if template_name in template_hit:
+    >>>             final_target_template_for_modeling[template_hit] = target_final_seq
+    >>>
+    >>> for template_pdb, target_seq in final_target_template_for_modeling.items():
+    >>>     output_model_name = 'protein_comp_modeling/{}_{}.pdb'.format(target_seq_path.split('.')[0], '_'.join(template_pdb.split('_')[0:2]))
+    >>>     join_apo_dir_path = os.path.join(template_pdb_path, template_pdb)
+    >>>     pose = pyrosetta.pose_from_file(join_apo_dir_path)
+    >>>     assert(pose.size() == len(target_seq))
+    >>>     scorefxn = pyrosetta.get_fa_scorefxn()
+    >>>     for i in range(len(target_seq)):
+    >>>         seqpos = i + 1
+    >>>         name1 = target_seq[i]
+    >>>         if (name1 == "-"):
+    >>>             continue
+    >>>         pyrosetta.rosetta.protocols.toolbox.pose_manipulation.repack_this_residue(seqpos, pose, scorefxn, True, name1)
+    >>>     pose.dump_pdb(output_model_name)
+    """

kinoml/workflows/kinase_ligand_modeling.py

@@ -0,0 +1,158 @@
+"""
+Workflows for kinase modeling, based on Karanicolas' protocols
+
+Given a ligand and a target kinase, this workflow will:
+
+Part A: preparation
+
+    1. Align the ligand against an in-house database of active kinase ligands
+       (resulting in a total of 100 conformers).
+    2. Align the kinase sequence against an in-house database of kinase sequences
+
+Part B1: build best models
+
+    3. Use Rosetta to build 10 most promising models, as suggested by the alignment
+       similarity.
+    4. Best model will be combined with conformers described in step 1 (100 total
+       complexes)
+    5. Minimize them.
+
+Part B2: build remaining models
+
+    6. Pick 10 best ligands from step 1, and all 10 models from step 10. Build the
+       resulting 100 models.
+    7. Minimize them.
+
+Part C: Results
+
+    8. Report results from B1 and B2
+
+
+Notes
+-----
+
+Originally developed at https://github.com/karanicolaslab/kinmodel.
+
+"""
+
+from argparse import ArgumentParser
+import os
+import numpy as np
+from openforcefield.topology import Molecule
+
+
+def parse_cli():
+    p = ArgumentParser()
+    p.add_argument("ligand")
+    p.add_argument("protein")
+    p.add_argument("--nconformers", default=100, help="Number of ligand conformers that will be generated")
+    p.add_argument("--nmodels", default=10, help="Number of target models that will be built")
+    p.add_argument("--ligand_library", default=None)
+    p.add_argument("--sequence_library", default=None)
+
+    return p.parse_args()
+
+
+def _load_ligand(ligand):
+    """
+    Load a ligand using OpenForceField toolkit.
+
+    Parameters
+    ----------
+    ligand : str
+        This can be a path to a file (supported by OFF) or a SMILES string
+
+    Returns
+    -------
+    openforcefield.topology.Molecule
+    """
+    if os.path.isfile(ligand):
+        return Molecule.from_file(ligand)
+    return Molecule.from_smiles(ligand)
+
+
+def _superpose_ligand(ligand, dataset):
+    """
+    Structural superposition of a query molecule ``ligand`` against an existing
+    dataset of small compounds
+
+    Parameters
+    ----------
+    ligand : openforcefield.topology.Molecule
+        Query ligand that will be superposed. Conformers must have been generated
+        previously.
+    dataset : str or list of openforcefield.topology.Molecule
+        Path to a SDF file containing the target database
+
+    Returns
+    -------
+    superposed_conformers : array, shape=len(dataset)*len(ligand.conformers)*ligand.n_atoms*3
+    scores : array, shape=len(dataset)*len(ligand.conformers)
+    """
+    if isinstance(dataset, str):
+        dataset = _load_dataset(dataset)
+    for molecule in dataset:
+        # align(molecule, ligand)
+        pass
+
+
+def _load_dataset(path):
+    """
+    Load a molecule dataset into a list of OFF Molecules
+
+    Parameters
+    ----------
+    path : str or list of str
+        Path to the dataset file (SDF, MOL) or multiple files.
+        Supported formats are SDF, MOL, PDB
+
+    Returns
+    -------
+    list of openforcefield.topology.Molecule
+    """
+    if isinstance(path, (list, tuple)):
+        pass
+    pass
+
+def _align_sequence():
+    pass
+
+
+def main():
+    # Temporarily add imports here to have a broader picture of the code
+
+    args = parse_cli()
+
+    # Part A1 - align ligands
+    ligand_molecule = _load_ligand(args.ligand)
+    ligand_molecule.generate_conformers(n_conformers=args.nconformers)
+
+    superposed_xyz, superposition_scores = _superpose_ligand(ligand_molecule, args.ligand_library)
+
+    # TODO: Get best conformers: see https://stackoverflow.com/a/38884051
+    best_conformers_indices = np.argpartition(superposition_scores, -args.nconformers)[-args.nconformers:]
+    best_conformers = superposed_xyz[best_conformers_indices]  # FIXME: This does not work as expected
+
+    # Part A2 - align protein sequence and build protein models
+    protein_molecule = _load_protein(args.protein)
+    protein_sequence = protein_molecule.sequence
+    aligned_sequences, alignment_scores = _align_sequence(protein_sequence, args.sequence_library)
+    best_alignment_indices = np.argpartition(alignment_scores, -args.nmodels)[-args.nmodels:]
+    candidate_sequences = aligned_sequences[best_alignment_indices]  # FIXME
+    built_models = build_with_rosetta(template, sequences, n_results=10)  # TODO: parallel call!
+
+    # Part B1 - build protein-ligand complexes and minimize
+    complexes = []
+    for protein, ligand in product(built_models, best_conformers):
+        complexes.append(concatenate(protein, ligand))
+
+    # TODO: parallel call
+    minimized_complexes = [minimize(c) for c in complexes]
+
+    # Part B2 - repeat for remaining models
+
+    # Part C - report
+    report(...)
+
+
+