From 4f5c317d84c15d0b39c4885bdc98793fd96a8d11 Mon Sep 17 00:00:00 2001 From: Jaime RGP Date: Mon, 16 Sep 2019 12:19:57 +0200 Subject: [PATCH 1/4] Add draft script to start PR --- kinoml/workflows/kinase_ligand_modeling.py | 128 +++++++++++++++++++++ 1 file changed, 128 insertions(+) create mode 100644 kinoml/workflows/kinase_ligand_modeling.py diff --git a/kinoml/workflows/kinase_ligand_modeling.py b/kinoml/workflows/kinase_ligand_modeling.py new file mode 100644 index 00000000..c092e007 --- /dev/null +++ b/kinoml/workflows/kinase_ligand_modeling.py @@ -0,0 +1,128 @@ +""" +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. + +""" + +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 + 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) + """ + + +def main(): + 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(...) + + + From b4904617ee10b9ed7a58a4ca2af957a6d33c8fdf Mon Sep 17 00:00:00 2001 From: Jaime RGP Date: Thu, 17 Oct 2019 07:01:27 +0200 Subject: [PATCH 2/4] Stash some temporary changes --- kinoml/workflows/kinase_ligand_modeling.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/kinoml/workflows/kinase_ligand_modeling.py b/kinoml/workflows/kinase_ligand_modeling.py index c092e007..481eadb1 100644 --- a/kinoml/workflows/kinase_ligand_modeling.py +++ b/kinoml/workflows/kinase_ligand_modeling.py @@ -88,9 +88,16 @@ def _superpose_ligand(ligand, dataset): superposed_conformers : array, shape=len(dataset)*len(ligand.conformers)*ligand.n_atoms*3 scores : array, shape=len(dataset)*len(ligand.conformers) """ + 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 From 68eb4d816da5d04de88a505ee02d5da29eaf074a Mon Sep 17 00:00:00 2001 From: Jaime RGP Date: Fri, 18 Oct 2019 16:09:07 +0200 Subject: [PATCH 3/4] add stubs in docstrings for further refactor --- kinoml/modeling/complexes.py | 18 +++++++ kinoml/modeling/ligands.py | 64 +++++++++++++++++++++++++ kinoml/modeling/receptors.py | 93 ++++++++++++++++++++++++++++++++++++ 3 files changed, 175 insertions(+) create mode 100644 kinoml/modeling/complexes.py create mode 100644 kinoml/modeling/ligands.py create mode 100644 kinoml/modeling/receptors.py diff --git a/kinoml/modeling/complexes.py b/kinoml/modeling/complexes.py new file mode 100644 index 00000000..4968346b --- /dev/null +++ b/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)) + """ \ No newline at end of file diff --git a/kinoml/modeling/ligands.py b/kinoml/modeling/ligands.py new file mode 100644 index 00000000..54c07e97 --- /dev/null +++ b/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) diff --git a/kinoml/modeling/receptors.py b/kinoml/modeling/receptors.py new file mode 100644 index 00000000..10ae09c8 --- /dev/null +++ b/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) + """ From b9172284ec643aa92504ef4023c7499860f1d2fa Mon Sep 17 00:00:00 2001 From: Jaime RGP Date: Wed, 23 Oct 2019 16:36:34 +0200 Subject: [PATCH 4/4] Add more stubs --- kinoml/workflows/kinase_ligand_modeling.py | 27 ++++++++++++++++++++-- 1 file changed, 25 insertions(+), 2 deletions(-) diff --git a/kinoml/workflows/kinase_ligand_modeling.py b/kinoml/workflows/kinase_ligand_modeling.py index 481eadb1..7b650162 100644 --- a/kinoml/workflows/kinase_ligand_modeling.py +++ b/kinoml/workflows/kinase_ligand_modeling.py @@ -35,6 +35,7 @@ """ +from argparse import ArgumentParser import os import numpy as np from openforcefield.topology import Molecule @@ -80,7 +81,7 @@ def _superpose_ligand(ligand, dataset): ligand : openforcefield.topology.Molecule Query ligand that will be superposed. Conformers must have been generated previously. - dataset : str + dataset : str or list of openforcefield.topology.Molecule Path to a SDF file containing the target database Returns @@ -88,8 +89,30 @@ def _superpose_ligand(ligand, dataset): superposed_conformers : array, shape=len(dataset)*len(ligand.conformers)*ligand.n_atoms*3 scores : array, shape=len(dataset)*len(ligand.conformers) """ - pass + 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