How to use lightdock - 10 common examples

#!/usr/bin/env python

"""Calculates few statistics about the result of the simulation"""

import sys
import os

from lightdock.util.logger import LoggingManager
from matplotlib.backends.backend_agg import FigureCanvasAgg as FigureCanvas
from matplotlib.figure import Figure

log = LoggingManager.get_logger('stats')


def usage():
    """Displays usage parameters and exists with error"""
    log.error("Wrong arguments")
    raise SystemExit("usage: %s number_of_steps number_of_glowworms" % (sys.argv[0]))


def parse_command_line():
    # Arguments parsing
    if len(sys.argv) != 3:
        usage()
    try:
        num_steps = int(sys.argv[1])
    except Exception, e:
        log.error(str(e))

neighbors = []
    vision_range = []
    scoring = []

    data_file = open(lightdock_output)
    lines = [line.rstrip(os.linesep) for line in data_file.readlines()]
    data_file.close()

    counter = 0
    for line in lines:
        if line[0] == '(':
            counter += 1
            last = line.index(')')
            coord = line[1:last].split(',')
            translations.append([float(coord[0]), float(coord[1]), float(coord[2])])
            rotations.append(Quaternion(float(coord[3]), float(coord[4]), float(coord[5]), float(coord[6])))
            values = line[last + 1:].split()
            luciferin.append(float(values[0]))
            neighbors.append(int(values[1]))
            vision_range.append(float(values[2]))
            scoring.append(float(values[3]))

    log.info("Read %s coordinate lines" % counter)
    return translations, rotations, luciferin, neighbors, vision_range, scoring

ligand_ids = []
    rec_extents = []
    lig_extents = []

    data_file = open(lightdock_output)
    lines = data_file.readlines()
    data_file.close()

    counter = 0
    for line in lines:
        if line[0] == '(':
            counter += 1
            last = line.index(')')
            coord = line[1:last].split(',')
            translations.append([float(coord[0]), float(coord[1]), float(coord[2])])
            rotations.append(Quaternion(float(coord[3]), float(coord[4]), float(coord[5]), float(coord[6])))
            if len(coord) > 7:
                rec_extents.append(np.array([float(x) for x in coord[7:7+num_anm_rec]]))
                lig_extents.append(np.array([float(x) for x in coord[-num_anm_lig:]]))
            raw_data = line[last+1:].split()
            receptor_id = int(raw_data[0])
            ligand_id = int(raw_data[1])
            receptor_ids.append(receptor_id)
            ligand_ids.append(ligand_id)
    log.info("Read %s coordinate lines" % counter)
    return translations, rotations, receptor_ids, ligand_ids, rec_extents, lig_extents

oda_file_name = DEFAULT_LIGHTDOCK_PREFIX % (molecule.structure_file_names[0] + '.oda')
        if os.path.exists(oda_file_name):
            log.info('ODA %s file found' % oda_file_name)
            with open(oda_file_name) as oda_input:
                oda_values = []
                lines = oda_input.readlines()[2:]
                for line in lines:
                    try:
                        fields = line.split()
                        oda_values.append(float(fields[3]))
                    except ValueError:
                        pass
        return oda_values


class SIPPER(ScoringFunction):
    def __init__(self, weight=1.0):
        super(SIPPER, self).__init__(weight)
        self.energy = sipper_energy

    def __call__(self, receptor, receptor_coordinates, ligand, ligand_coordinates):
        """Computes the pyDock scoring energy using receptor and ligand which are
        instances of DockingModel
        """
        energy = csipper.calculate_sipper(receptor_coordinates, ligand_coordinates, self.energy,
                                          receptor.indexes, ligand.indexes,
                                          receptor.atoms_per_residue, ligand.atoms_per_residue,
                                          len(receptor.atoms_per_residue), len(ligand.atoms_per_residue),
                                          receptor.oda, ligand.oda)
        return energy * self.weight

data_file = open(data_file_name)
        data = data_file.readlines()
        data_file.close()
        potentials = [[None for i in range(22)] for j in range(22)]  # @UnusedVariable
        try:
            for x in range(22):
                for y in range(22):
                    potentials[x][y] = float(data[x + 1].strip().split()[y + 1])

        except Exception, e:
            raise PotentialsParsingError('Error parsing %s file. Details: %s' % (data_file_name, str(e)))

        return potentials


class TOBIAdapter(ModelAdapter):
    """Adapts a given Complex to a DockingModel object suitable for this
    TOBI scoring function.
    """

    def _get_docking_model(self, molecule, restraints):
        """Builds a suitable docking model for this scoring function"""
        residues = [residue for chain in molecule.chains for residue in chain.residues]
        tobi_residues = []
        list_of_coordinates = []
        parsed_restraints = {}
        for structure in range(molecule.num_structures):
            coordinates = []
            for residue in residues:
                try:
                    residue_index = TOBIPotential.recognized_residues.index(residue.name)

parser.add_argument("receptor", help="PDB receptor")
    parser.add_argument("ligand", help="PDB ligand")
    script_args = parser.parse_args()
    return script_args


if __name__ == "__main__":
    args = parse_command_line()

    try:
        scoring_function_module = "lightdock.scoring.%s.driver" % args.scoring_function
        module = importlib.import_module(scoring_function_module)
    except ImportError:
        raise SystemExit("Scoring function not found or not available")

    atoms, residues, chains = parse_complex_from_file(args.receptor)
    receptor = Complex(chains, atoms, structure_file_name=args.receptor)
    atoms, residues, chains = parse_complex_from_file(args.ligand)
    ligand = Complex(chains, atoms, structure_file_name=args.ligand)

    CurrentScoringFunction = getattr(module, "DefinedScoringFunction")
    CurrentModelAdapter = getattr(module, "DefinedModelAdapter")

    adapter = CurrentModelAdapter(receptor, ligand)
    scoring_function = CurrentScoringFunction()

    energy = scoring_function(adapter.receptor_model,  adapter.receptor_model.coordinates[0],
                              adapter.ligand_model, adapter.ligand_model.coordinates[0])
    print args.scoring_function, ': ', energy

args = parse_command_line()

        atoms_to_ignore = []
        if args.noxt:
            atoms_to_ignore.append('OXT')

        structures = []
        file_names = []
        file_name, file_extension = os.path.splitext(args.structure)
        if file_extension == DEFAULT_LIST_EXTENSION:
            file_names.extend(get_pdb_files(args.structure))
        else:
            file_names.append(args.structure)
        for file_name in file_names:
            log.info("Reading %s PDB file..." % file_name)
            atoms, residues, chains = parse_complex_from_file(file_name, atoms_to_ignore)
            structures.append({'atoms': atoms, 'residues': residues, 'chains': chains, 'file_name': file_name})
            log.info("%s atoms, %s residues read." % (len(atoms), len(residues)))

        molecule = Complex.from_structures(structures)
        try:
            ellipsoid = MinimumVolumeEllipsoid(molecule.atom_coordinates[0].coordinates)
        except MinimumVolumeEllipsoidError, e:
            log.error("Impossible to calculate minimum volume ellipsoid. Reason: %s" % str(e))
            raise SystemExit("%s finished with error" % script_name)

        output_file_name = molecule.structure_file_names[0] + DEFAULT_REFERENCE_POINTS_EXTENSION
        with open(output_file_name, 'w') as output:
            for point in ellipsoid.poles:
                output.write(get_point_respresentation(point) + os.linesep)
            output.write(get_point_respresentation(ellipsoid.center) + os.linesep)
        log.info('Points written to %s' % output_file_name)

script_args = parser.parse_args()
    return script_args


if __name__ == "__main__":
    args = parse_command_line()

    try:
        scoring_function_module = "lightdock.scoring.%s.driver" % args.scoring_function
        module = importlib.import_module(scoring_function_module)
    except ImportError:
        raise SystemExit("Scoring function not found or not available")

    atoms, residues, chains = parse_complex_from_file(args.receptor)
    receptor = Complex(chains, atoms, structure_file_name=args.receptor)
    atoms, residues, chains = parse_complex_from_file(args.ligand)
    ligand = Complex(chains, atoms, structure_file_name=args.ligand)

    CurrentScoringFunction = getattr(module, "DefinedScoringFunction")
    CurrentModelAdapter = getattr(module, "DefinedModelAdapter")

    adapter = CurrentModelAdapter(receptor, ligand)
    scoring_function = CurrentScoringFunction()

    energy = scoring_function(adapter.receptor_model,  adapter.receptor_model.coordinates[0],
                              adapter.ligand_model, adapter.ligand_model.coordinates[0])
    print args.scoring_function, ': ', energy

atoms_to_ignore.append('OXT')

        structures = []
        file_names = []
        file_name, file_extension = os.path.splitext(args.structure)
        if file_extension == DEFAULT_LIST_EXTENSION:
            file_names.extend(get_pdb_files(args.structure))
        else:
            file_names.append(args.structure)
        for file_name in file_names:
            log.info("Reading %s PDB file..." % file_name)
            atoms, residues, chains = parse_complex_from_file(file_name, atoms_to_ignore)
            structures.append({'atoms': atoms, 'residues': residues, 'chains': chains, 'file_name': file_name})
            log.info("%s atoms, %s residues read." % (len(atoms), len(residues)))

        molecule = Complex.from_structures(structures)
        try:
            ellipsoid = MinimumVolumeEllipsoid(molecule.atom_coordinates[0].coordinates)
        except MinimumVolumeEllipsoidError, e:
            log.error("Impossible to calculate minimum volume ellipsoid. Reason: %s" % str(e))
            raise SystemExit("%s finished with error" % script_name)

        output_file_name = molecule.structure_file_names[0] + DEFAULT_REFERENCE_POINTS_EXTENSION
        with open(output_file_name, 'w') as output:
            for point in ellipsoid.poles:
                output.write(get_point_respresentation(point) + os.linesep)
            output.write(get_point_respresentation(ellipsoid.center) + os.linesep)
        log.info('Points written to %s' % output_file_name)

        points = [point for point in ellipsoid.poles]
        points.append(ellipsoid.center)
        pdb_file_name = output_file_name + '.pdb'

n_modes = int(sys.argv[2])
        factor = float(sys.argv[3])
    except:
        usage()
        raise SystemExit("Wrong command line")

    protein = parsePDB(pdb_structure)
    ca_atoms = protein.select('name CA')
    protein_anm = ANM('protein ca')
    protein_anm.buildHessian(ca_atoms)
    protein_anm.calcModes(n_modes=n_modes)
    print 'Normal modes calculated'

    atoms, residues, chains = parse_complex_from_file(pdb_structure)
    lightdock_structures = [{'atoms': atoms, 'residues': residues, 'chains': chains, 'file_name': pdb_structure}]
    lightdock_structure = Complex.from_structures(lightdock_structures)
    print 'Structure read by lightdock'

    num_atoms_prody = len(protein.protein)
    num_atoms_lightdock = len(atoms)

    if num_atoms_prody != num_atoms_lightdock:
        raise SystemExit("Number of atoms is different")

    protein_anm_ext, protein_all = extendModel(protein_anm, ca_atoms, protein, norm=True)

    modes = []
    for i in range(n_modes):
        nm = protein_anm_ext.getEigvecs()[:, i].reshape((num_atoms_lightdock, 3))
        modes.append(nm)

    coordinates = lightdock_structure.atom_coordinates[0].coordinates