Commit ef26e227 authored by Evan N. Feinberg's avatar Evan N. Feinberg
Browse files

Merge pull request #57 from evanfeinberg/master

Added nnscore featurization and nosetest
parents b7490970 242ff9b9
Loading
Loading
Loading
Loading
+93 −0
Original line number Diff line number Diff line
"""
Contains methods for generating a pdbbind dataset mapping
  complexes (protein + ligand) to experimental binding measurement.
"""
import pickle
import os
import pandas as pd
from rdkit import Chem
from glob import glob
import re


def extract_labels(pdbbind_label_file):
  """Extract labels from pdbbind label file."""
  assert os.path.isfile(pdbbind_label_file)
  labels = {}
  with open(pdbbind_label_file) as f:
    content = f.readlines()
    for line in content:
      if line[0] == "#":
        continue
      line = line.split()
      # lines in the label file have format
      # PDB-code Resolution Release-Year -logKd Kd reference ligand-name
      #print line[0], line[3]
      labels[line[0]] = line[3]
  return labels

def construct_df(pdb_stem_directory, pdbbind_label_file, pdbbind_df_pkl):
  """
  Takes as input a stem directory containing subdirectories with ligand
    and protein pdb/mol2 files, a pdbbind_label_file containing binding
    assay data for the co-crystallized ligand in each pdb file,
    and a pdbbind_df_pkl to which will be saved a pandas DataFrame
    where each row contains a pdb_id, smiles string, unique complex id,
    ligand pdb as a list of strings per line in file, protein pdb as a list
    of strings per line in file, ligand mol2 as a list of strings per line in
    mol2 file, and a "label" containing the experimental measurement.
  """
  labels = extract_labels(pdbbind_label_file)
  df_rows = []
  os.chdir(pdb_stem_directory)
  pdb_directories = [pdb.replace('/', '') for pdb in glob('*/')]

  for pdb_dir in pdb_directories:
    print "About to extract ligand and protein input files"
    pdb_id = os.path.basename(pdb_dir)
    ligand_pdb = None
    protein_pdb = None
    for f in os.listdir(pdb_dir):
      if re.search("_ligand_hyd.pdb$", f):
        ligand_pdb = f
      elif re.search("_protein_hyd.pdb$", f):
        protein_pdb = f
      elif re.search("_ligand.mol2$", f):
        ligand_mol2 = f

    print "Extracted Input Files:"
    print (ligand_pdb, protein_pdb, ligand_mol2)
    if not ligand_pdb or not protein_pdb or not ligand_mol2:
      raise ValueError("Required files not present for %s" % pdb_dir)
    ligand_pdb_path = os.path.join(pdb_dir, ligand_pdb)
    protein_pdb_path = os.path.join(pdb_dir, protein_pdb)
    ligand_mol2_path = os.path.join(pdb_dir, ligand_mol2)

    with open(protein_pdb_path, "rb") as f:
      protein_pdb_lines = f.readlines()

    with open(ligand_pdb_path, "rb") as f:
      ligand_pdb_lines = f.readlines()

    try:
      with open(ligand_mol2_path, "rb") as f:
        ligand_mol2_lines = f.readlines()
    except:
      ligand_mol2_lines = []

    print "About to compute ligand smiles string."
    ligand_mol = Chem.MolFromPDBFile(ligand_pdb_path)
    if ligand_mol is None:
      continue
    smiles = Chem.MolToSmiles(ligand_mol)
    complex_id = "%s%s" % (pdb_id, smiles)
    label = labels[pdb_id]
    df_rows.append([pdb_id, smiles, complex_id, protein_pdb_lines,
                    ligand_pdb_lines, ligand_mol2_lines, label])

  pdbbind_df = pd.DataFrame(df_rows, columns=('pdb_id', 'smiles', 'complex_id',
                                              'protein_pdb', 'ligand_pdb',
                                              'ligand_mol2', 'label'))

  with open(pdbbind_df_pkl, "wb") as f:
    pickle.dump(pdbbind_df, f)
+32.1 MiB

File added.

No diff preview for this file type.

+5 −2
Original line number Diff line number Diff line
@@ -128,11 +128,12 @@ class Model(object):
    for epoch in range(self.model_params["nb_epoch"]):
      print("Starting epoch %s" % str(epoch+1))
      for i, (X, y, w, _) in enumerate(dataset.itershards()):
        print("Training on batch-%s/epoch-%s" % (str(i+1), str(epoch+1)))
        print("Training on shard-%s/epoch-%s" % (str(i+1), str(epoch+1)))
        nb_sample = np.shape(X)[0]
        interval_points = np.linspace(
            0, nb_sample, np.ceil(float(nb_sample)/batch_size)+1, dtype=int)
        for j in range(len(interval_points)-1):
          print("Training on batch-%s/shard-%s/epoch-%s" % (str(j+1), str(i+1), str(epoch+1)))
          indices = range(interval_points[j], interval_points[j+1])
          X_batch = X[indices, :]
          y_batch = y[indices]
@@ -160,7 +161,9 @@ class Model(object):
      y_preds = []
      for j in range(len(interval_points)-1):
        indices = range(interval_points[j], interval_points[j+1])
        y_preds.append(self.predict_on_batch(X[indices, :]))
        y_pred_on_batch = self.predict_on_batch(X[indices, :]).reshape((len(indices),len(task_names)))
        y_preds.append(y_pred_on_batch)

      y_pred = np.concatenate(y_preds)
      y_pred = np.reshape(y_pred, np.shape(y))

+30 −17
Original line number Diff line number Diff line
@@ -45,8 +45,14 @@ def add_featurize_group(featurize_cmd):
      "--threshold", type=float, default=None,
      help="If specified, will be used to binarize real-valued target-fields.")
  featurize_group.add_argument(
      "--parallel", type=float, default=None,
      help="Use multiprocessing will be used to parallelize featurization.")
      "--protein-pdb-field", type=str, default=None,
      help="Name of field holding protein pdb.")
  featurize_group.add_argument(
      "--ligand-pdb-field", type=str, default=None,
      help="Name of field holding ligand pdb.")
  featurize_group.add_argument(
      "--ligand-mol2-field", type=str, default=None,
      help="Name of field holding ligand mol2.")

def add_transforms_group(cmd):
  """Adds flags for data transforms."""
@@ -62,7 +68,7 @@ def add_transforms_group(cmd):
           "to mean no transforms are required.")
  transform_group.add_argument(
      "--feature-types", nargs="+", required=1,
      choices=["user-specified-features", "ECFP", "RDKIT-descriptors"],
      choices=["user-specified-features", "ECFP", "RDKIT-descriptors", "NNScore"],
      help="Featurizations of data to use.\n"
           "'features' denotes user-defined features.\n"
           "'fingerprints' denotes ECFP fingeprints.\n"
@@ -205,7 +211,8 @@ def create_model(args):
    featurize_inputs(
        feature_dir, data_dir, args.input_files, args.user_specified_features,
        args.tasks, args.smiles_field, args.split_field, args.id_field,
        args.threshold, args.parallel)
        args.threshold, args.protein_pdb_field,
        args.ligand_pdb_field, args.ligand_mol2_field)

  if args.generate_dataset:
    print("+++++++++++++++++++++++++++++++++")
@@ -273,9 +280,14 @@ def parse_args(input_args=None):
  add_model_command(subparsers)
  return parser.parse_args(input_args)

def shard_inputs(input_file):
  input_file_no_ext = os.path.splitext(input_file)


def featurize_inputs(feature_dir, data_dir, input_files,
                     user_specified_features, tasks, smiles_field,
                     split_field, id_field, threshold, parallel):
                     split_field, id_field, threshold, protein_pdb_field, 
                     ligand_pdb_field, ligand_mol2_field):

  """Allows for parallel data featurization."""
  featurize_input_partial = partial(featurize_input,
@@ -285,13 +297,11 @@ def featurize_inputs(feature_dir, data_dir, input_files,
                                    smiles_field=smiles_field,
                                    split_field=split_field,
                                    id_field=id_field,
                                    threshold=threshold)
                                    threshold=threshold,
                                    protein_pdb_field=protein_pdb_field,
                                    ligand_pdb_field=ligand_pdb_field,
                                    ligand_mol2_field=ligand_mol2_field)

  if parallel:
    pool = mp.Pool(int(mp.cpu_count()/2))
    pool.map(featurize_input_partial, input_files)
    pool.terminate()
  else:
  for input_file in input_files:
    featurize_input_partial(input_file)

@@ -302,18 +312,21 @@ def featurize_inputs(feature_dir, data_dir, input_files,
  FeaturizedSamples(samples_dir, dataset_files)

def featurize_input(input_file, feature_dir, user_specified_features, tasks,
                    smiles_field, split_field, id_field, threshold):
                    smiles_field, split_field, id_field, threshold, protein_pdb_field,
                     ligand_pdb_field, ligand_mol2_field):
  """Featurizes raw input data."""
  featurizer = DataFeaturizer(tasks=tasks,
                              smiles_field=smiles_field,
                              split_field=split_field,
                              id_field=id_field,
                              threshold=threshold,
                              protein_pdb_field=protein_pdb_field,
                              ligand_pdb_field=ligand_pdb_field,
                              ligand_mol2_field=ligand_mol2_field,
                              user_specified_features=user_specified_features,
                              verbose=True)
  out = os.path.join(
      feature_dir, "%s.joblib" %(os.path.splitext(os.path.basename(input_file))[0]))
  featurizer.featurize(input_file, FeaturizedSamples.feature_types, out)

  featurizer.featurize(input_file, FeaturizedSamples.feature_types, feature_dir)

def train_test_split(input_transforms, output_transforms,
                     feature_types, splittype, data_dir):
+2 −1
Original line number Diff line number Diff line
@@ -314,7 +314,8 @@ def _df_to_numpy(df, feature_types):
  sorted_ids = df["mol_id"]

  # Set missing data to have weight zero
  missing = (y == "")
  missing = (y.astype(object) == "")

  y[missing] = 0.
  w[missing] = 0.

Loading