[med-svn] [plip] 03/05: New upstream version 1.3.2+dfsg

Alex Mestiashvili malex-guest at moszumanska.debian.org
Mon Sep 19 15:01:49 UTC 2016


This is an automated email from the git hooks/post-receive script.

malex-guest pushed a commit to branch master
in repository plip.

commit 2b99b24c066d096c9fad811d556f3dcb4b3db4e9
Author: Alexandre Mestiashvili <alex at biotec.tu-dresden.de>
Date:   Mon Sep 19 16:55:44 2016 +0200

    New upstream version 1.3.2+dfsg
---
 CHANGES.txt                            |   7 +
 DOCUMENTATION.md                       | 365 +++++++++++++++++++++++++++------
 README.md                              | 105 +++++++---
 plip/modules/config.py                 |   4 +-
 plip/modules/detection.py              |   5 +-
 plip/modules/preparation.py            | 111 +++++++---
 plip/modules/pymolplip.py              |  10 +-
 plip/modules/report.py                 |  87 +++++++-
 plip/modules/supplemental.py           |  11 +
 plip/modules/visualize.py              |  43 +++-
 plip/modules/webservices.py            |  67 ++++++
 plip/plipcmd                           | 121 ++---------
 plip/test/test_basic_functions.py      |  90 +++++++-
 plip/test/test_literature_validated.py |   4 +-
 plip/test/test_remote_services.py      |  46 +++++
 setup.py                               |   2 +-
 16 files changed, 828 insertions(+), 250 deletions(-)

diff --git a/CHANGES.txt b/CHANGES.txt
index 431fb7b..921ce28 100644
--- a/CHANGES.txt
+++ b/CHANGES.txt
@@ -1,6 +1,13 @@
 Changelog
 ---------
 
+### 1.3.2
+* __Processing of protein-peptide interactions via --peptides__
+* option to keep modified residues as ligands (--keepmod)
+* Improved code for reports
+* Smart ordering of ligand in composite compounds
+* Fixes handling and visualization of DNA/RNA
+
 ### 1.3.1
 * __Support for amino acids as ligands__
 * __Plugin-ready for PyMOL and Chimera__
diff --git a/DOCUMENTATION.md b/DOCUMENTATION.md
index 7099037..947cc0d 100644
--- a/DOCUMENTATION.md
+++ b/DOCUMENTATION.md
@@ -1,122 +1,345 @@
-===========
-PLIP v1.2.2
-===========
+# Protein-Ligand Interaction Profiler (PLIP)
 
+## What is PLIP?
 The Protein-Ligand Interaction Profiler (PLIP) is a tool to analyze and visualize protein-ligand interactions in PDB files.
 
+## Why should I use PLIP?
 
-Features
-========
-* Detection of eight different types of noncovalent interactions, including metal complexes
-* Works for complexes of protein with small molecules, ions, polymers, or DNA/RNA (and all combinations)
+#### Comprehensive Detection of Interaction
+* Eight types of noncovalent interactions
+* Interaction between proteins and
+  * small molecules
+  * ions
+  * polymers
+  * DNA / RNA
+* Rich additiona information on binding, e.g. unpaired functional groups
+
+#### Everything Is Automated
+
+* Direct download of PDB structures from PDB server
 * Automatic detection and grouping of relevant ligands in a PDB file
-* Rich additional information on binding, e.g. unpaired functional groups
-* Direct download of PDB structures from wwPDB server if valid PDB ID is given
-* Processing of custom PDB files containing protein-ligand complexes (e.g. from docking)
 * No need for special preparation of a PDB file, works out of the box
+* Automatic fixiing of errors in PDB files
+
+#### Flexible Usage
+* Processing of custom PDB files containing protein-ligand complexes (e.g. from docking)
 * Atom-level interaction reports in TXT and XML formats for easy parsing
 * Generation of PyMOL session files (.pse), enabling easy preparation of images for publications and talks
 * Rendering of 3D interaction diagram for each ligand and its interactions with the protein
 
-Quickstart
-==========
-Run the PLIP python script with the option `-h` to list all available command line options.
-To analyze a protein-ligand complex from a Protein Data Bank entry -- e.g. 1vsn --, run
-    `plipcmd -i 1vsn`.
-To analyze a PDB file from your workstation, run
-    `plipcmd -f path_to_pdbfile.pdb`.
-The output format(s) can be chosen freely, ranging from ...
-* XML report files (`-x`, highest level of detail)
-* Text report files (`-t`, medium level of detail)
+## Get PLIP running
+
+#### 1. Install dependencies
+The following instructions work on a unix machine.
+If you are using Windows or OSX, the process may differ.
+Please refer to the webpages of the corresponding tools to get help for installation.
+
+To install all dependencies for PLIP, either just install PyMOL and then use the `pip` installation routine or install all tools and clone from GitHub.
+The current version of PLIP depends on
+* OpenBabel >=2.3.2
+* PyMOL 1.7.x with Python bindings
+* Imagemagick >=6.9.x (optional)
+
+and should be executed with Python 2.7.x.
+
+Example command for Ubuntu using apt-get
+```bash
+sudo apt-get install pymol openbabel python-openbabel imagemagick
+```
+
+
+#### 2. Get PLIP on your machine
+
+There are several options to get PLIP on your machine.
+If you have the Python package manager `pip` installed, you can install PLIP in a terminal via the following command:
+
+```bash
+pip install plip
+```
+If you want to clone from Github, open a new system terminal and clone the repository using
+```bash
+git clone https://github.com/ssalentin/plip.git ~/pliptool
+```
+You can use any other user directory, but we will use the one given above for the documentation.
+
+#### 3. Simplify access to PLIP
+
+The package manager `pip` will create symlinks for you so you can call PLIP using `plip` in the command line. If you cloned from Github, use the following command to do the same:
+
+```bash
+alias plip='~/pliptool/plip/plipcmd'
+```
+
+## Analyze a PDB structure with PLIP
+
+### Single Structures
+Having PLIP installed, you can run
+```bash
+plip -h
+```
+to list all available command line options.
+We will go through all important settings in this section.
+A typical application of PLIP involves the analysis of protein-ligand interactions in a structure from the Protein Data Bank (PDB).
+PLIP can automatically fetch the entry from the PDB server when a valid PDB ID is provided.
+```bash
+plip -i 1vsn -v
+```
+The command above will fetch the PDB entry 1vsn from the server, analyze all interactions and print out the results (verbose mode).
+No output files are produced at this point.
+
+The same can be done for local PDB files.
+
+```bash
+wget http://files.rcsb.org/download/1EVE.pdb
+plip -f 1EVE.pdb -v
+```
+
+The output formats can be added in any combination, currently including:
+* XML report files (`-x`, best for automatic processing)
+* Text report files (`-t`, human-readable)
 * PyMOL session files (`-y`)
 * Ray-traced images (`-p`)
-* Verbose output on command line (`-v`)
 
-Threshold settings
-==================
-All geometric thresholds used for detection of interactions in PLIP are stored in the `config.py` module and can be
-changed there permanently if desired. Another possibility is to change single parameters via command line options when
-running plip. The naming of the variables is identical to those in the config file, but all lowercase. Specify the
+```bash
+plip -i 1osn -pyx
+```
+
+The command above will fetch the PDB entry 1osn, analyze all interactions and produce one XML result file (`-x`) as well as rendered images (`-p`) and PyMOL session files (`-y`) for each binding site.
+
+### Batch Mode
+PLIP can process multiple structures at once, either from PDB or local files.
+To activate batch mode, just provide a list of PDB IDs or local file names, e.g.:
+
+```bash
+plip -i 1vsn 1osn 2reg -vx
+```
+
+PLIP will create subdirectories for each given structure in the output folder.
+If in PDB ID mode (`i`), the folder structure will be nested and based on the two middle characters of the PDB ID.
+The structure 1vsn in batch processing will have its output files in <outputfolder>/vs/1vsn .
+
+### Detection of protein-peptide interactions
+For the detection of ligands, PLIP relies on the separation of ATOM and HETATM entries in the PDB file.
+The latter are searched for suitable ligands when running in normal mode.
+Peptide ligands, however, are usually deposited as ATOM entries in a separate chain.
+PLIP can not detect these entities automatically.
+To switch into protein-peptide interaction mode, start PLIP with the option `--peptide`, followed by the peptide chain of interest, e.g.:
+
+```bash
+plip -i 5hi4 --peptides I -vx
+```
+
+## Changing detection thresholds
+The PLIP algorithm uses a rule-based detection to report non-covalent interaction between proteins and their partners.
+The current settings are based on literature values and have been refined based on extensive testing with independent cases from mainly crystallography journals, covering a broad range of structure resolutions.
+For most users, it is not recommended to change the standard settings.
+However, there may be cases where changing detection thresholds is advisable (e.g. sets with multiple very low resolution structures)
+PLIP allows you to change the setting permanently or for one run.
+
+#### Permanent change
+PLIP settings are stored in the file `modules/config.py`, which is loaded in each run.
+One section in the file includes assignments of standard values to all distance and angle thresholds, together with short descriptions.
+Changing the values in this file means they will be changed permanently for PLIP.
+Furthermore, filter settings for ligands (metal ions in complexes, unsupported ligands, warnings for possible artifacts) can be changed in this file as well.
+
+#### Temporary change
+
+Thresholds can be changed for each run using command-line parameters.
+The naming of the variables is identical to those in the config file, but all lowercase.
+Specify the
 threshold you want to change and the new value, e.g. to change HYDROPH_DIST_MAX to 5 Angstrom, run PLIP using
-    `python plip-cmd.py -i 1vsn --hydroph_dist_max 5`
+```bash
+plip -i 1vsn --hydroph_dist_max 5
+```
 All distance thresholds can be increased to up to 10 Angstrom. Thresholds for angles can be set between 0 and 180 degree.
 If two interdependent thresholds have conflicting values, PLIP will show an error message.
 
-Web Service
-===========
+## Further Options
+PLIP offers further command line options which enables you to switch advanced settings.
+* Set number `n` of maximum threads used for parallel processing (`--maxthreads <n>`)
+* Do not automatically combine covalently bound ligands (`--breakcomposite`)
+* Do not discard alternate locations (`--altlocation`)
+* Set debug mode (`--debug`)
+* Turn off automatic fixing of errors in PDB files (`--nofix`)
+* Keep modified residues as ligands (`--keepmod`)
+
+## Web Service
 A web service for analysis of protein-ligand complexes using PLIP is available at
-http://projects.biotec.tu-dresden.de/plip-web
+http://plip.biotec.tu-dresden.de/
 The web site offers advanced functions to search for specific entries from PDB and lists the interaction results in the browser.
+Additionally, the service used the BioLiP database to annotate biologically relevant ligands.
+The option to change threshold, ligand filtering, and batch processing is only available in the command line tool and with the Python modules.
+
+## Algorithm
+PLIP uses a rule-based system for detection of non-covalent interactions between protein residues and ligands.
+Information on chemical groups able to participate in a specific interaction (e.g. requirements for hydrogen bond donors) and interaction geometry (e.g. distance and angle thresholds) from literature are used to detect characteristics of non-covalent interactions between contacting atoms of protein and ligands.
+For each binding site, the algorithm searches first for atoms or atom groups in the protein and ligand which could possibly be partner in specific interactions.
+In the second step, geometric rules are applied to match groups in protein and ligand forming an interaction.
+
+#### Detection and filtering of ligands
+Previous to the detection step for the interactions, PLIP extracts all ligands contained in the structure.
+Modified amino acids are identified and excluded using MODRES entries of the PDB files.
+Additionally, it uses the BioLiP list of possible artifacts to remove ligands which are in this list and appear 15 times or more in a structure.
+Just a few compounds are currently excluded, being listed in the PLIP config file.
+
+#### Preparation of structures
+Polar hydrogens are added to the structure and alternative conformations/models/positions removed.
+Missing chains are assigned to ligands and non-standard ligand names (with special characters) altered to LIG.
+
+#### Detection of possible interacting groups
+
+##### Binding site atoms
+The binding site distance cutoff is determined by adding up BS_DIST_MAX to the maximum extent to the ligand (maximum distance of a ligand atom to ligand centroid).
+All protein atoms within this distance cutoff to any binding site atoms are counted as belonging to the binding site.
+
+##### Hydrophobic Atoms
+An atom is classified as hydrophobic if it is a carbon and has only carbon or hydrogen atoms as neighbours.
+
+##### Aromatic Rings
+OpenBabel is used to identify rings (SSSR perception) and their aromaticity.
+In cases where no aromaticity is reported by OpenBabel, the ring is checked for planarity.
+To this end, the normals of each atom in the ring to its neighbors is calculated.
+The angle between each pair of normals has to be less than AROMATIC_PLANARITY.
+If this holds true, the ring is also considered as aromatic.
+
+##### Hydrogen Bond Donors and Acceptors
+OpenBabel is used to identify hydrogen bond donor and acceptor atoms.
+Halogen atoms are excluded from this group and treated separately (see below).
+
+##### Charged Groups
+The detection of charged groups is only exhaustive for the binding site, not the ligands.
+For proteins, positive charges are attributed to the side chain nitrogens of Arginine, Histidine and Lysine.
+Negative charged are assigned to the carboxyl groups in Aspartic Acid and Glutamic Acid.
+In ligands, positive charges are assigned to quaterny ammonium groups, tertiary amines (assuming the nitrogen could pick up a hydrogen and thus get charged), sulfonium and guanidine groups.
+Negative charges are reported for phosphate, sulfonate, sulfonic acid and carboxylate.
+
+##### Halogen bonds donors and acceptors
+Assuming that halogen atoms are not present in proteins (unless they are artificially modified), halogen bond donors are searched for only in ligands.
+All fluorine, chlorine, bromide or iodine atoms connected to a carbon atom qualify as donors.
+Halogen bond acceptors in proteins are all carbon, phosphor or sulphur atoms connected to oxygen, phosphor, nitrogen or sulfur.
 
+##### Water
+Water atoms are assigned to a ligand-binding site complex if their oxygen atoms are within a certain cutoff to the ligand.
+The cutoff is determined by adding up BS_DIST_MAX to the maximum extent to the ligand (maximum distance of a ligand atom to ligand centroid).
+This means the farthest distance of a ligand to a water atom is BS_DIST_MAX.
+
+#### Detection of interactions
+For an overview on geometric cutoffs used for the prediction of interactions, please refer to the config file.
+Note that the threshold can not be changed for jobs running on PLIP Web.
+The command line tool (sourcecode available for download), however, allows changing all listed parameters permanently or for single runs.
+
+##### Hydrophobic Interactions
+As hydrophobic interactions result from entropic changes rather than attractive forces between atoms, there are no clear geometries of hydrophobic association.
+The observed attraction between hydrophobic atoms decays exponentionally with the distance between them.
+
+A generous cutoff was chosen, identifying a prime set of hydrophobic interactions between all pairs of hydrophobic atoms within a distance of HYDROPH_DIST_MAX.
+Since the number of hydrophobic interactions with such an one-step approach can easily surpass all other interaction types combined, it may strongly influence subsequent evaluation or applications as interaction fingerprinting.
+To overcome this problem, the number of hydrophobic interactions is reduced in several steps.
+First, hydrophobic interactions between rings interacting via π-stacking are removed.
+This is done because stacking already involves hydrophobic interactions.
 
-Documentation
-=============
+Second, two clustering steps are applied.
+If a ligand atom interacts with several binding site atoms in the same residue, only the interaction with the closest distance is kept.
+Subsequently, the set of hydrophobic interactions is checked from the opposite perspective: if a protein atom interacts with several neighboring ligand atoms, just the interaction with the closest distance is kept.
+Together, these reduction steps help to report only the most representative hydrophobic interactions.
 
-Algorithm
----------
-PLIP uses a rule-based algorithm for the detection of non-covalent interactions. For details on the algorithm, visit
-http://projects.biotec.tu-dresden.de/plip-web
+##### Hydrogen Bonds
+A hydrogen bond between a hydrogen bond donor and acceptor is reported if several geometric requirements are fulfilled.
+The distance has to be less than HBOND_DIST_MAX and the angle at the donor group (D-H...A) above HBOND_DON_ANGLE_MIN.
 
-Exit codes
-----------
-1 : Unspecified Error
-2 : Empty PDB file as input
-3 : Invalid PDB ID (wrong format)
-4 : PDB file can't be read by OpenBabel (due to invalid input files)
-5 : PDB ID is valid, but wwPDB offers no file in PDB format for download.
+Since salt bridges involve purely electrostatic interactions as well as hydrogen bonds, it is not meaningful to report both interaction types between the same groups.
+Thus, hydrogen bonds between atoms are removed if they belong to groups that already form a salt bridge to that atom.
 
-Legend for PyMOL visualization
-------------------------------
-All colors given as RGB values.
-<Description> - <RGB> - <PyMOL color> - <Representation>
+As a general rule, a hydrogen bond donor can take part in only one hydrogen bond, while acceptor atoms can be partners in multiple hydrogen bonds (e.g. bifurcated hydrogen bonds).
+For multiple possible hydrogen bonds from one donor, only the contact with the donor angle closer to 180 ° is kept.
 
-__Structural Elements__
+##### π-Stacking
+π-Stacking for two aromatic rings is reported whenever their centers are within a distance of PISTACK_DIST_MAX, the angle deviates no more than PISTACK_ANG_DEV from the optimal angle of 90 ° for T-stacking or 180 ° for P-stacking.
 
+Additionally, each ring center is projected onto the opposite ring plane.
+The distance between the other ring center and the projected point (i.e. the offset) has to be less than PISTACK_OFFSET_MAX.
+This value corresponds approximately to the radius of benzene + 0.6 Å.
 
-Protein - [43, 131, 186] - myblue (custom) - sticks
+##### π-Cation Interactions
+π-Cation interactions are reported for each pairing of a positive charge and an aromatic ring if the distance between the charge center and the aromatic ring center is less than PICATION_DIST_MAX.
+In the case of a putative π-cation interaction with a tertiary amine of the ligand, an additional angle criterion is applied (see documentation in the source code).
 
-Ligand - [253, 174, 97] - myorange (custom) - sticks
+##### Salt Bridges
+Whenever two centers of opposite charges come within a distance of SALTBRIDGE_DIST_MAX, a salt bridge is reported.
+In contrast to hydrogen bonds, there are no additonal geometric restrictions.
 
-Water - [191, 191, 255] - lightblue - nb_spheres
+##### Water bridges
+While residues can be bridged by more than one water molecule, for the prediction in this script the only case considered is one water molecule bridging ligand and protein atoms via hydrogen bonding.
 
-Charge Center - [255, 255, 0] - yellow - spheres
+The water molecule has to be positioned between hydrogen bond donor/acceptor pairs of ligand and protein with distances of the water oxygen within WATER_BRIDGE_MINDIST and WATER_BRIDGE_MAXDIST to the corresponding polar atoms of the donor or acceptor groups.
+If a constellation with a water atom fulfils these requirements, two angles are checked.
+The angle ω between the acceptor atom, the water oxygen and donor hydrogen has to be within WATER_BRIDGE_OMEGA_MIN and WATER_BRIDGE_OMEGA_MAX.
+Additionally, the angle θ between the water oxygen, the donor hydrogen and the donor atom has to be larger than WATER_BRIDGE_THETA_MIN.
 
-Aromatic Ring Center - [230, 230, 230] -  grey90 - spheres
+Similar to standard hydrogen bonds, a water molecule is only allowed to participate as donor in two hydrogen bonds (two hydrogen atoms as donors).
+In the case of more than two possible hydrogen bonds for a water molecule as donor, only the two contacts with a water angle closest to 110 ° are kept
 
-Ions - [255, 255, 128] - hotpink - spheres
+##### Halogen Bonds
+Halogen bonds are reported for each pairing of halogen bond acceptor and donor group having a distance of less than HALOGEN_DIST_MAX and angles at the donor and acceptor group of HALOGEN_DON_ANGLE and HALOGEN_ACC_ANGLE with a deviation of no more than HALOGEN_ANG_DEV
 
+##### Metal Complexes
+For metal complexes, PLIP considers metal ions from a set of more than 50 species (see PLIP config for more details).
+Possible interacting groups in the protein are sidechains of cystein (S), histidine (N), asparagine, glutamic acid, serin, threonin, and tyrosin (all O), as well as all main chain oxygens.
 
-__Interactions__
+In ligands, following groups are considered for metal complexation: alcohols, phenolates, carboxylates, phosphoryls, thiolates, imidazoles, pyrroles, and the iron-sulfur cluster as a special constellation.
+For one metal ions, all groups with a maximum distance of METAL_DIST_MAX to the ligand are considered for the complex.
 
+After assigning all target groups to one metal ions, the resulting set of angles of the complex is compared with known sets of angles from common coordination geometries (linear [2], trigonal planar [3], trigonal pyramidal [3], tetrahedral [4], square planar [4], trigonal bipyramidal [5], square pyramidal [5], and octahedral [6]).
+The best fit with the least difference in observed targets is chosen as an estimated geometry and targets superfluous to the constellation are removed.
 
-Hydrophobic Interaction - [128, 128, 128] - grey50 - dashed Line
+## Additional Information
 
-Hydrogen Bond - [0, 0, 255] - blue - solid Line
+### Exit error codes
 
-Water Bridges - [191, 191, 255] - lightblue - solid Line
+| Exit code  | Description |
+| ------------- | ------------- |
+| **1**  | Unspecified Error  |
+| **2**  | Empty PDB file as input |
+| **3**  | Invalid PDB ID (wrong format)  |
+| **4**  | PDB file can't be read by OpenBabel  |
+| **5**  | Valid PDB ID, but no PDB file on wwwPDB  |
 
-pi-Stacking (parallel) - [0, 255, 0] - green - dashed Line
+### Legend for PyMOL visualization
 
-pi-Stacking (perpendicular) - [140, 179, 102] - smudge - dashed Line
+#### Structural Elements
 
-pi-Cation Interaction - [255, 128, 0] - orange - dashed Line
+| Description  | RGB | PyMOL color | Representation |
+| ------------ | --- | ------------| ---------------|
+| Protein  | [43, 131, 186] | myblue (custom) | sticks |
+| Ligand  | [253, 174, 97] | myorange (custom) | sticks |
+| Water  | [191, 191, 255]  | lightblue | nb_spheres |
+| Charge Center | [255, 255, 0] | yellow | spheres |
+| Aromatic Ring Center  | [230, 230, 230] | grey90 | spheres |
+| Ions | [250, 255, 128] | hotpink | spheres |
 
-Halogen Bond - [64, 255, 191] - greencyan - solid Line
 
-Salt Bridge - [255, 255, 0] - yellow - dashed Line
+#### Interactions
 
-Metal Complexation - [140, 64, 153] - violetpurple - dashed Line
+| Description  | RGB | PyMOL color | Representation |
+| ------------ | --- | ------------| -------------- |
+| Hydrophobic Interaction  | [128, 128, 128] | grey50 | dashed |
+| Hydrogen Bond  | [0, 0, 255] | blue | solid line |
+| Water Bridges  | [191, 191, 255]  | lightblue | solid line |
+| pi-Stacking (parallel) | [0, 255, 0] | green | dashed line |
+| pi-Stacking (perpendicular)  | [140, 179, 102] | smudge | dashed line |
+| pi-Cation Interaction | [255, 128, 0] | orange | dashed line |
+| Halogen Bond | [54, 255, 191] | greencyan | solid line |
+| Salt Bridge | [255, 255, 0] | yellow | dashed line |
+| Metal Complex | [140, 64, 153] | violetpurple | dashed line |
 
 
-XML Report Documentation
-------------------------
 
+### XML Report Documentation
 
 **report**
 
 Contains all binding site information
 
-
 **plipversion**
 
 Version of PLIP used for generating the output file
@@ -260,3 +483,11 @@ Coordinates of interacting ligand atom or interaction center in ligand
 **protcoo**
 
 Coordinates of interacting protein atom or interaction center in ligand
+
+## Testing
+To run the tests, run the following set of commands.
+
+```bash
+cd ~/pliptool/plip/test
+python -m unittest test*
+```
diff --git a/README.md b/README.md
index fb8cb63..df1d8cf 100644
--- a/README.md
+++ b/README.md
@@ -1,43 +1,94 @@
-PLIP
-====
+# PLIP
 
-Protein-Ligand Interaction Profiler (PLIP) - Analyze non-covalent protein-ligand interactions in 3D structures
+**Protein-Ligand Interaction Profiler (PLIP)**
 
-Versions
---------
-The latest commits may contain newer, but untested features. The last version marked as release is always stable.
 
-Installation
-------------
+ Analyze non-covalent protein-ligand interactions in 3D structures
 
-Previous to the installation of PLIP, make sure you have the following tools and libraries installed:
-* Python 2.x
+## How to use PLIP
+
+#### 1. Clone the repository
+
+Open a new system terminal and clone this repository using
+```bash
+git clone https://github.com/ssalentin/plip.git ~/pliptool
+```
+#### 2. Run PLIP
+
+##### As a command line tool
+
+Run the `plipcmd` script inside the PLIP folder to detect, report, and visualize interactions. The following example creates a PYMOL visualization for the interactions
+between the inhibitor NFT and its target protein in the PDB structure 1VSN.
+
+```bash
+alias plip='~/pliptool/plip/plipcmd'
+mkdir /tmp/1vsn && cd /tmp/1vsn
+plip -i 1vsn -yv
+pymol 1VSN_NFT_A_283.pse
+```
+
+##### As a python library
+
+In your terminal, add the PLIP repository to your PYTHONPATH variable.
+For our example, we also download a PDB file for testing.
+```bash
+export PYTHONPATH=~/pliptool/plip:${PYTHONPATH}
+cd /tmp && wget http://files.rcsb.org/download/1EVE.pdb
+python
+```
+In python, import the PLIP modules, load a PDB structure and run the analysis.
+This small example shows how to print all numbers of residues involved in pi-stacking:
+
+```python
+from plip.modules.preparation import PDBComplex
+
+my_mol = PDBComplex()
+my_mol.load_pdb('/tmp/1EVE.pdb') # Load the PDB file into PLIP class
+print my_mol # Shows name of structure and ligand binding sites
+my_bsid = 'E20:A:2001' # Unique binding site identifier (HetID:Chain:Position)
+my_mol.analyze()
+my_interactions = my_mol.interaction_sets[my_bsid] # Contains all interaction data
+
+# Now print numbers of all residues taking part in pi-stacking
+print [pistack.resnr for pistack in s.pistacking] # Prints [84, 129]
+```
+
+#### 3. View and process the results
+PLIP offers various output formats, ranging from renderes images and PyMOL session files to human-readable text files and XML files.
+By default, all files are deposited in the working directory unless and output path is provided.
+For a full documentation of running options and output formats, please refear to the documentation.
+
+## Versions and Branches
+For production environments, you should use the latest versioned commit from the stable branch.
+Newer commits from the stable and development branch may contains new but untested and not documented features.
+
+## Requirements
+Previous to using PLIP, make sure you have the following tools and libraries installed:
+* Python 2.7.x
 * OpenBabel >=2.3.2
 * PyMOL 1.7.x with Python bindings
-* pip (Python package managaer)
-* Imagemagick (optional, needed for automatic cropping of images)
-
-To install PLIP, simply run the following command using a terminal
-> python setup.py install
+* Imagemagick >=6.9.x (optional)
 
-Code Contributions
-------------------
+## Contributions
 Sebastian Salentin sebastian.salentin (at) biotec.tu-dresden.de
 
-Joachim Haupt joachim.haupt (at) biotec.tu-dresden.de
+Joachim Haupt joachim.haupt (at) biotec.tu-dresden.de | https://github.com/vjhaupt
+
+Melissa F. Adasme Mora melissa.adasme (at) biotec.tu-dresden.de
 
-Melissa F. Adasme Mora melissa.adasme (at) biotec.tu-dresden.de (Testing and Validation)
+## PLIP Web Server
+Visit our PLIP Web Server on http://plip.biotec.tu-dresden.de/
 
-PLIP Web Server
----------------
-Visit our PLIP Web Server on http://projects.biotec.tu-dresden.de/plip-web
+## Contact Me
+Do you have feature requests, found a bug or want to use `PLIP` in your project?
+Write me an email to sebastian.salentin (at) biotec.tu-dresden.de
 
-Contact Me
-----------
-Questions or comments about `PLIP`? Write me an email to sebastian.salentin (at) biotec.tu-dresden.de
+## License Information
+PLIP is published under the Apache License. For more information, please read the LICENSE.txt file.
+Using PLIP in your commercial or non-commercial project is generally possible when giving a proper reference to this project and the publication in NAR.
+If you are unsure about usage options, don't hesitate to contact me.
 
-Citation Information
---------------------
+## Citation Information
 If you are using PLIP in your work, please cite
 > Salentin,S. et al. PLIP: fully automated protein-ligand interaction profiler.
 > Nucl. Acids Res. (1 July 2015) 43 (W1): W443-W447. doi: 10.1093/nar/gkv315
diff --git a/plip/modules/config.py b/plip/modules/config.py
index 7f78ad4..3f1a8ea 100644
--- a/plip/modules/config.py
+++ b/plip/modules/config.py
@@ -23,13 +23,14 @@ XML = False
 TXT = False
 PICS = False
 PYMOL = False
-VISJSON = False  # JSON File of PyMOL helper class (for web service)
 OUTPATH = './'
 BASEPATH = './'
 BREAKCOMPOSITE = False  # Break up composite ligands with covalent bonds
 ALTLOC = False  # Consider alternate locations
 PLUGIN_MODE = False  # Special mode for PLIP in Plugins (e.g. PyMOL)
 NOFIX = False  # Turn off fixing of errors in PDB files
+PEPTIDES = []  # Definition which chains should be considered as peptide ligands
+KEEPMOD = False
 
 # Configuration file for Protein-Ligand Interaction Profiler (PLIP)
 # Set thresholds for detection of interactions
@@ -79,6 +80,7 @@ other = ['MO', 'RE', 'HO']
 UNSUPPORTED = anions + other
 
 # BioLiP list of suspicious ligands from http://zhanglab.ccmb.med.umich.edu/BioLiP/ligand_list (2014-07-10)
+# Add ligands here to get warnings for possible artifacts.
 biolip_list = ['ACE', 'HEX', 'TMA', 'SOH', 'P25', 'CCN', 'PR', 'PTN', 'NO3', 'TCN', 'BU1', 'BCN', 'CB3', 'HCS', 'NBN',
                'SO2', 'MO6', 'MOH', 'CAC', 'MLT', 'KR', '6PH', 'MOS', 'UNL', 'MO3', 'SR', 'CD3', 'PB', 'ACM', 'LUT',
                'PMS', 'OF3', 'SCN', 'DHB', 'E4N', '13P', '3PG', 'CYC', 'NC', 'BEN', 'NAO', 'PHQ', 'EPE', 'BME', 'TB',
diff --git a/plip/modules/detection.py b/plip/modules/detection.py
index 5aa64f9..2087552 100644
--- a/plip/modules/detection.py
+++ b/plip/modules/detection.py
@@ -372,8 +372,9 @@ def metal_complexation(metals, metal_binding_lig, metal_binding_bs):
                     final_geom, final_coo, = next_total.geometry, next_total.coordination
                     rms, excluded = next_total.rms, next_total.excluded
                     break
-                elif i == len(all_total) - 1:
-                    final_geom, final_coo, rms, excluded = "NA", "NA", 0.0, []
+                elif i == len(all_total) - 2:
+                    final_geom, final_coo, rms, excluded = "NA", "NA", float('nan'), []
+                    break
 
         # Record all contact pairing, excluding those with targets superfluous for chosen geometry
         only_water = set([x[0].location for x in contact_pairs]) == {'water'}
diff --git a/plip/modules/preparation.py b/plip/modules/preparation.py
index d083c9b..a8b12be 100644
--- a/plip/modules/preparation.py
+++ b/plip/modules/preparation.py
@@ -182,41 +182,72 @@ class LigandFinder:
         self.ligands = self.getligs()
         self.excluded = sorted(list(self.lignames_all.difference(set(self.lignames_kept))))
 
+    def getpeptides(self, chain):
+        """If peptide ligand chains are defined via the command line options,
+        try to extract the underlying ligand formed by all residues in the
+        given chain without water
+        """
+        data = namedtuple('ligand', 'mol hetid chain position water members longname type atomorder can_to_pdb')
+        all_from_chain = [o for o in pybel.ob.OBResidueIter(self.proteincomplex.OBMol) if o.GetChain() == chain] # All residues from chain
+        if len(all_from_chain) == 0:
+            return None
+        else:
+            #TODO Can be done shorter with a split
+            water = [o for o in all_from_chain if o.GetResidueProperty(9)]
+            non_water = water = [o for o in all_from_chain if not o.GetResidueProperty(9)]
+            ligand = self.extract_ligand(non_water)
+            return ligand
+
+
     def getligs(self):
         """Get all ligands from a PDB file and prepare them for analysis.
         Returns all non-empty ligands.
         """
-        ligands = []
 
-        # Filter for ligands using lists
-        ligand_residues, self.lignames_all, self.water = self.filter_for_ligands()
+        if config.PEPTIDES == []:
+            # Extract small molecule ligands (default)
+            ligands = []
+
+            # Filter for ligands using lists
+            ligand_residues, self.lignames_all, self.water = self.filter_for_ligands()
 
-        all_res_dict = {(a.GetName(), a.GetChain(), a.GetNum()): a for a in ligand_residues}
-        self.lignames_kept = list(set([a.GetName() for a in ligand_residues]))
+            all_res_dict = {(a.GetName(), a.GetChain(), a.GetNum()): a for a in ligand_residues}
+            self.lignames_kept = list(set([a.GetName() for a in ligand_residues]))
+
+            if not config.BREAKCOMPOSITE:
+                #  Update register of covalent links with those between DNA/RNA subunits
+                self.covalent += nucleotide_linkage(all_res_dict)
+                #  Find fragment linked by covalent bonds
+                res_kmers = self.identify_kmers(all_res_dict)
+            else:
+                res_kmers = [[a, ] for a in ligand_residues]
+            for kmer in res_kmers:  # iterate over all ligands and extract molecules + information
+                ligands.append(self.extract_ligand(kmer))
 
-        if not config.BREAKCOMPOSITE:
-            #  Update register of covalent links with those between DNA/RNA subunits
-            self.covalent += nucleotide_linkage(all_res_dict)
-            #  Find fragment linked by covalent bonds
-            res_kmers = self.identify_kmers(all_res_dict)
         else:
-            res_kmers = [[a, ] for a in ligand_residues]
-        for kmer in res_kmers:  # iterate over all ligands and extract molecules + information
-            ligands.append(self.extract_ligand(kmer))
+            # Extract peptides from given chains
+            self.water = [o for o in pybel.ob.OBResidueIter(self.proteincomplex.OBMol) if o.GetResidueProperty(9)]
+            peptide_ligands = [self.getpeptides(chain) for chain in config.PEPTIDES]
+            ligands = [p for p in peptide_ligands if p is not None ]
+            self.covalent, self.lignames_kept, self.lignames_all = [], [], set()
+
         return [lig for lig in ligands if len(lig.mol.atoms) != 0]
 
     def extract_ligand(self, kmer):
         """Extract the ligand by copying atoms and bonds and assign all information necessary for later steps."""
         data = namedtuple('ligand', 'mol hetid chain position water members longname type atomorder can_to_pdb')
         members = [(res.GetName(), res.GetChain(), int32_to_negative(res.GetNum())) for res in kmer]
-        members = sorted(members, key=lambda x: (x[1], x[2]))
+        members = sort_members_by_importance(members)
         rname, rchain, rnum = members[0]
         write_message("Finalizing extraction for ligand %s:%s:%s\n" % (rname, rchain, rnum), mtype='debug')
         names = [x[0] for x in members]
         longname = '-'.join([x[0] for x in members])
 
-        # Classify a ligand by its HETID(s)
-        ligtype = classify_by_name(names)
+        if config.PEPTIDES == []:
+            # Classify a ligand by its HETID(s)
+            ligtype = classify_by_name(names)
+        else:
+            ligtype = 'PEPTIDE'
 
         hetatoms = set()
         for obresidue in kmer:
@@ -275,29 +306,37 @@ class LigandFinder:
         return ligand
 
     def is_het_residue(self, obres):
-        """Given an OBResidue, determines if the residue is indeed a ligand
-        in the PDB file (any atoms has to be a HETATM entry)"""
-        het_atoms = []
-        for atm in pybel.ob.OBResidueAtomIter(obres):
-                het_atoms.append(obres.IsHetAtom(atm))
-        if True in het_atoms:
+        """Given an OBResidue, determines if the residue is indeed a possible ligand
+        in the PDB file"""
+        if not obres.GetResidueProperty(0):
+            # If the residue is NOT amino (0)
+            # It can be amino_nucleo, coenzme, ion, nucleo, protein, purine, pyrimidine, solvent
+            # In these cases, it is a ligand candidate
             return True
         else:
-            return False
+            # Here, the residue is classified as amino
+            # Amino acids can still be ligands, so we check for HETATM entries
+            # Only residues with at least one HETATM entry are processed as ligands
+            het_atoms = []
+            for atm in pybel.ob.OBResidueAtomIter(obres):
+                    het_atoms.append(obres.IsHetAtom(atm))
+            if True in het_atoms:
+                return True
+        return False
 
 
     def filter_for_ligands(self):
         """Given an OpenBabel Molecule, get all ligands, their names, and water"""
 
-        #candidates1 = [o for o in pybel.ob.OBResidueIter(self.proteincomplex.OBMol) if not (o.GetResidueProperty(9)
-        #                                                                                or o.GetResidueProperty(0))]
-
         candidates1 = [o for o in pybel.ob.OBResidueIter(self.proteincomplex.OBMol) if not o.GetResidueProperty(9) and self.is_het_residue(o)]
         all_lignames = set([a.GetName() for a in candidates1])
 
         water = [o for o in pybel.ob.OBResidueIter(self.proteincomplex.OBMol) if o.GetResidueProperty(9)]
         # Filter out non-ligands
-        candidates2 = [a for a in candidates1 if is_lig(a.GetName()) and a.GetName() not in self.modresidues]
+        if not config.KEEPMOD:  # Keep modified residues as ligands
+            candidates2 = [a for a in candidates1 if is_lig(a.GetName()) and a.GetName() not in self.modresidues]
+        else:
+            candidates2 = [a for a in candidates1 if is_lig(a.GetName())]
         write_message("%i ligand(s) after first filtering step.\n" % len(candidates2), mtype='debug')
 
         ############################################
@@ -1152,6 +1191,10 @@ class PDBComplex:
         self.Mapper = Mapper()
         self.ligands = []
 
+    def __str__(self):
+        formatted_lig_names = [":".join([x.hetid, x.chain, str(x.position)]) for x in self.ligands]
+        return "Protein structure %s with ligands:\n" % (self.pymol_name) + "\n".join([lig for lig in formatted_lig_names])
+
     def load_pdb(self, pdbpath, as_string=False):
         """Loads a pdb file with protein AND ligand(s), separates and prepares them.
         If specified 'as_string', the input is a PDB string instead of a path."""
@@ -1224,6 +1267,11 @@ class PDBComplex:
         else:
             write_message("Structure contains no ligands.\n\n")
 
+    def analyze(self):
+        """Triggers analysis of all complexes in structure"""
+        for ligand in self.ligands:
+            self.characterize_complex(ligand)
+
     def characterize_complex(self, ligand):
         """Handles all basic functions for characterizing the interactions for one ligand"""
 
@@ -1277,9 +1325,14 @@ class PDBComplex:
         return [obres.GetIdx() for obres in resis if self.res_belongs_to_bs(obres, cutoff, ligcentroid)]
 
     def res_belongs_to_bs(self, res, cutoff, ligcentroid):
-        """Check for each residue if its centroid is within a certain distance to the ligand centroid."""
+        """Check for each residue if its centroid is within a certain distance to the ligand centroid.
+        Additionally checks if a residue belongs to a chain restricted by the user (e.g. by defining a peptide chain)"""
         rescentroid = centroid([(atm.x(), atm.y(), atm.z()) for atm in pybel.ob.OBResidueAtomIter(res)])
-        return True if euclidean3d(rescentroid, ligcentroid) < cutoff else False
+        # Check geometry
+        near_enough = True if euclidean3d(rescentroid, ligcentroid) < cutoff else False
+        # Check chain membership
+        restricted_chain = True if res.GetChain() in config.PEPTIDES else False
+        return (near_enough and not restricted_chain)
 
     def get_atom(self, idx):
         return self.atoms[idx]
diff --git a/plip/modules/pymolplip.py b/plip/modules/pymolplip.py
index 4012dcc..14a9914 100644
--- a/plip/modules/pymolplip.py
+++ b/plip/modules/pymolplip.py
@@ -16,6 +16,7 @@ See the License for the specific language governing permissions and
 limitations under the License.
 """
 from pymol import cmd
+from time import sleep
 import sys
 import os
 import subprocess
@@ -26,7 +27,8 @@ class PyMOLVisualizer:
         if not plcomplex is None:
             self.plcomplex = plcomplex
             self.protname = plcomplex.pdbid  # Name of protein with binding site
-            self.ligname = plcomplex.hetid  # Name of ligand
+            self.hetid = plcomplex.hetid
+            self.ligname = "Ligand_" + self.hetid  # Name of ligand
             self.metal_ids = plcomplex.metal_ids
 
     def set_initial_representations(self):
@@ -328,9 +330,11 @@ class PyMOLVisualizer:
                 cmd.zoom(self.ligname, 3)
         cmd.origin(self.ligname)
 
-    def save_session(self, outfolder):
+    def save_session(self, outfolder, override=None):
         """Saves a PyMOL session file."""
-        filename = '%s_%s' % (self.protname.upper(), "_".join([self.ligname, self.plcomplex.chain, self.plcomplex.position]))
+        filename = '%s_%s' % (self.protname.upper(), "_".join([self.hetid, self.plcomplex.chain, self.plcomplex.position]))
+        if override is not None:
+            filename = override
         cmd.save("/".join([outfolder, "%s.pse" % filename]))
 
     def png_workaround(self, filepath, width=1200, height=800):
diff --git a/plip/modules/report.py b/plip/modules/report.py
index a5fd480..d8cbeb3 100644
--- a/plip/modules/report.py
+++ b/plip/modules/report.py
@@ -20,13 +20,96 @@ limitations under the License.
 # Python Standard Library
 import time
 from operator import itemgetter
+import sys
+
 
 # External libraries
 import lxml.etree as et
 
+__version__ = '1.3.2'
+
+class StructureReport:
+    """Creates reports (xml or txt) for one structure/"""
+    def __init__(self, mol):
+        self.mol = mol
+        self.excluded = self.mol.excluded
+        self.xmlreport = self.construct_xml_tree()
+        self.txtreport = self.construct_txt_file()
+        self.get_bindingsite_data()
+        self.outpath = mol.output_path
+
+    def construct_xml_tree(self):
+        """Construct the basic XML tree"""
+        report = et.Element('report')
+        plipversion = et.SubElement(report, 'plipversion')
+        plipversion.text = __version__
+        date_of_creation = et.SubElement(report, 'date_of_creation')
+        date_of_creation.text = time.strftime("%Y/%m/%d")
+        citation_information = et.SubElement(report, 'citation_information')
+        citation_information.text = "Salentin,S. et al. PLIP: fully automated protein-ligand interaction profiler. " \
+           "Nucl. Acids Res. (1 July 2015) 43 (W1): W443-W447. doi: 10.1093/nar/gkv315"
+        pdbid = et.SubElement(report, 'pdbid')
+        pdbid.text = self.mol.pymol_name.upper()
+        filetype = et.SubElement(report, 'filetype')
+        filetype.text = self.mol.filetype.upper()
+        pdbfile = et.SubElement(report, 'pdbfile')
+        pdbfile.text = self.mol.sourcefiles['pdbcomplex']
+        pdbfixes = et.SubElement(report, 'pdbfixes')
+        pdbfixes.text = str(self.mol.information['pdbfixes'])
+        filename = et.SubElement(report, 'filename')
+        filename.text = str(self.mol.sourcefiles['filename'])
+        exligs = et.SubElement(report, 'excluded_ligands')
+        for i, exlig in enumerate(self.excluded):
+            e = et.SubElement(exligs, 'excluded_ligand', id=str(i + 1))
+            e.text = exlig
+        return report
+
+    def construct_txt_file(self):
+        """Construct the header of the txt file"""
+        textlines = ['Prediction of noncovalent interactions for PDB structure %s' % self.mol.pymol_name.upper(), ]
+        textlines.append("=" * len(textlines[0]))
+        textlines.append('Created on %s using PLIP v%s\n' % (time.strftime("%Y/%m/%d"), __version__))
+        textlines.append('If you are using PLIP in your work, please cite:')
+        textlines.append('Salentin,S. et al. PLIP: fully automated protein-ligand interaction profiler.')
+        textlines.append('Nucl. Acids Res. (1 July 2015) 43 (W1): W443-W447. doi: 10.1093/nar/gkv315\n')
+        if len(self.excluded) != 0:
+            textlines.append('Excluded molecules as ligands: %s\n' % ','.join([lig for lig in self.excluded]))
+        return textlines
+
+    def get_bindingsite_data(self):
+        """Get the additional data for the binding sites"""
+        for i, site in enumerate(sorted(self.mol.interaction_sets)):
+            s = self.mol.interaction_sets[site]
+            bindingsite = BindingSiteReport(s).generate_xml()
+            bindingsite.set('id', str(i + 1))
+            bindingsite.set('has_interactions', 'False')
+            self.xmlreport.insert(i + 1, bindingsite)
+            for itype in BindingSiteReport(s).generate_txt():
+                self.txtreport.append(itype)
+            if not s.no_interactions:
+                bindingsite.set('has_interactions', 'True')
+            else:
+                self.txtreport.append('No interactions detected.')
+            sys.stdout = sys.__stdout__  # Change back to original stdout, gets changed when PyMOL has been used before
+
+    def write_xml(self, as_string=False):
+        """Write the XML report"""
+        if not as_string:
+            et.ElementTree(self.xmlreport).write('%s/report.xml' % self.outpath, pretty_print=True, xml_declaration=True)
+        else:
+            print et.tostring(self.xmlreport, pretty_print=True)
+
+    def write_txt(self, as_string=False):
+        """Write the TXT report"""
+        if not as_string:
+            with open('%s/report.txt' % self.outpath, 'w') as f:
+                [f.write(textline + '\n') for textline in self.txtreport]
+        else:
+            print '\n'.join(self.txtreport)
+
 
-class TextOutput:
-    """Gather report data and generate reports for one binding site in different formats"""
+class BindingSiteReport:
+    """Gather report data and generate reports for one binding site in different formats."""
     def __init__(self, plcomplex):
 
         ################
diff --git a/plip/modules/supplemental.py b/plip/modules/supplemental.py
index 27398af..f49c538 100644
--- a/plip/modules/supplemental.py
+++ b/plip/modules/supplemental.py
@@ -318,6 +318,17 @@ def classify_by_name(names):
                     ligtype += '+ION'
     return ligtype
 
+def sort_members_by_importance(members):
+    """Sort the members of a composite ligand according to two criteria:
+    1. Split up in main and ion group. Ion groups are located behind the main group.
+    2. Within each group, sort by chain and position."""
+    main = [x for x in members if x[0] not in config.METAL_IONS]
+    ion = [x for x in members if x[0] in config.METAL_IONS]
+    sorted_main = sorted(main, key=lambda x: (x[1], x[2]))
+    sorted_main = sorted(main, key=lambda x: (x[1], x[2]))
+    sorted_ion = sorted(ion, key=lambda x: (x[1], x[2]))
+    return sorted_main + sorted_ion
+
 
 def get_isomorphisms(reference, lig):
     """Get all isomorphisms of the ligand."""
diff --git a/plip/modules/visualize.py b/plip/modules/visualize.py
index 5a72279..0de5ae7 100644
--- a/plip/modules/visualize.py
+++ b/plip/modules/visualize.py
@@ -18,7 +18,7 @@ limitations under the License.
 
 
 # Own modules
-from supplemental import initialize_pymol, start_pymol, write_message, colorlog
+from supplemental import initialize_pymol, start_pymol, write_message, colorlog, sysexit
 import config
 from pymolplip import PyMOLVisualizer
 from plipremote import VisualizerData
@@ -29,6 +29,7 @@ import sys
 
 # Special imports
 from pymol import cmd
+import pymol
 
 
 def select_by_ids(selname, idlist, selection_exists=False, chunksize=20, restrict=None):
@@ -60,7 +61,12 @@ def visualize_in_pymol(plcomplex):
     pdbid = plcomplex.pdbid
     lig_members = plcomplex.lig_members
     chain = plcomplex.chain
-    ligname = plcomplex.hetid
+    if config.PEPTIDES != []:
+        vis.ligname = 'PeptideChain%s' % plcomplex.chain
+
+    ligname = vis.ligname
+    hetid = plcomplex.hetid
+
     metal_ids = plcomplex.metal_ids
     metal_ids_str = '+'.join([str(i) for i in metal_ids])
 
@@ -76,7 +82,12 @@ def visualize_in_pymol(plcomplex):
     write_message('Setting current_name to "%s" and pdbid to "%s\n"' % (current_name, pdbid), mtype='debug')
     cmd.set_name(current_name, pdbid)
     cmd.hide('everything', 'all')
-    cmd.select(ligname, 'resn %s and chain %s and resi %s*' % (ligname, chain, plcomplex.position))
+    if config.PEPTIDES != []:
+        cmd.select(ligname, 'chain %s and not resn HOH' % plcomplex.chain)
+    else:
+        cmd.select(ligname, 'resn %s and chain %s and resi %s*' % (hetid, chain, plcomplex.position))
+    write_message("Selecting ligand for PDBID %s and ligand name %s with: " % (pdbid, ligname), mtype='debug')
+    write_message('resn %s and chain %s and resi %s*' % (hetid, chain, plcomplex.position), mtype='debug')
 
     # Visualize and color metal ions if there are any
     if not len(metal_ids) == 0:
@@ -84,9 +95,11 @@ def visualize_in_pymol(plcomplex):
         cmd.show('spheres', 'id %s and %s' % (metal_ids_str, pdbid))
 
     # Additionally, select all members of composite ligands
-    for member in lig_members:
-        resid, chain, resnr = member[0], member[1], str(member[2])
-        cmd.select(ligname, '%s or (resn %s and chain %s and resi %s)' % (ligname, resid, chain, resnr))
+    if len(lig_members) > 1:
+        for member in lig_members:
+           resid, chain, resnr = member[0], member[1], str(member[2])
+           cmd.select(ligname, '%s or (resn %s and chain %s and resi %s)' % (ligname, resid, chain, resnr))
+
     cmd.show('sticks', ligname)
     cmd.color('myblue')
     cmd.color('myorange', ligname)
@@ -97,6 +110,7 @@ def visualize_in_pymol(plcomplex):
         cmd.set('sphere_scale', 0.3, ligname)
     cmd.deselect()
 
+
     vis.make_initial_selections()
 
     vis.show_hydrophobic()  # Hydrophobic Contacts
@@ -116,8 +130,15 @@ def visualize_in_pymol(plcomplex):
     vis.selections_cleanup()
     vis.selections_group()
     vis.additional_cleanup()
-    if config.PYMOL:
-        vis.save_session(config.OUTPATH)
-    if config.PICS:
-        filename = '%s_%s' % (pdbid.upper(), "_".join([ligname, plcomplex.chain, plcomplex.position]))
-        vis.save_picture(config.OUTPATH, filename)
+    if config.PEPTIDES == []:
+        if config.PYMOL:
+            vis.save_session(config.OUTPATH)
+        if config.PICS:
+            filename = '%s_%s' % (pdbid.upper(), "_".join([hetid, plcomplex.chain, plcomplex.position]))
+            vis.save_picture(config.OUTPATH, filename)
+    else:
+        filename = "%s_PeptideChain%s" % (pdbid.upper(), plcomplex.chain)
+        if config.PYMOL:
+            vis.save_session(config.OUTPATH, override=filename)
+        if config.PICS:
+            vis.save_picture(config.OUTPATH, filename)
diff --git a/plip/modules/webservices.py b/plip/modules/webservices.py
new file mode 100644
index 0000000..c1f283f
--- /dev/null
+++ b/plip/modules/webservices.py
@@ -0,0 +1,67 @@
+"""
+Protein-Ligand Interaction Profiler - Analyze and visualize protein-ligand interactions in PDB files.
+webservices.py - Connect to various webservices to retrieve data
+Copyright 2014-2016 Sebastian Salentin
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+"""
+
+# Own modules
+from supplemental import write_message, sysexit
+
+# Python standard libary
+import urllib2
+
+# External libraries
+import lxml.etree as et
+
+
+
+def check_pdb_status(pdbid):
+    """Returns the status and up-to-date entry in the PDB for a given PDB ID"""
+    url = 'http://www.rcsb.org/pdb/rest/idStatus?structureId=%s' % pdbid
+    xmlf = urllib2.urlopen(url)
+    xml = et.parse(xmlf)
+    xmlf.close()
+    status = None
+    current_pdbid = pdbid
+    for df in xml.xpath('//record'):
+        status = df.attrib['status']  # Status of an entry can be either 'UNKWOWN', 'OBSOLETE', or 'CURRENT'
+        if status == 'OBSOLETE':
+            current_pdbid = df.attrib['replacedBy']  # Contains the up-to-date PDB ID for obsolete entries
+    return [status, current_pdbid.lower()]
+
+
+def fetch_pdb(pdbid):
+    """Get the newest entry from the RCSB server for the given PDB ID. Exits with '1' if PDB ID is invalid."""
+    pdbid = pdbid.lower()
+    write_message('\nChecking status of PDB ID %s ... ' % pdbid)
+    state, current_entry = check_pdb_status(pdbid)  # Get state and current PDB ID
+
+    if state == 'OBSOLETE':
+        write_message('entry is obsolete, getting %s instead.\n' % current_entry)
+    elif state == 'CURRENT':
+        write_message('entry is up to date.\n')
+    elif state == 'UNKNOWN':
+        sysexit(3, 'Invalid PDB ID (Entry does not exist on PDB server)\n')
+    write_message('Downloading file from PDB ... ')
+    pdburl = 'http://www.rcsb.org/pdb/files/%s.pdb' % current_entry  # Get URL for current entry
+    try:
+        pdbfile = urllib2.urlopen(pdburl).read()
+        # If no PDB file is available, a text is now shown with "We're sorry, but ..."
+        # Could previously be distinguished by an HTTP error
+        if 'sorry' in pdbfile:
+            sysexit(5, "No file in PDB format available from wwPDB for the given PDB ID.\n")
+    except urllib2.HTTPError:
+        sysexit(5, "No file in PDB format available from wwPDB for the given PDB ID.\n")
+    return [pdbfile, current_entry]
diff --git a/plip/plipcmd b/plip/plipcmd
index 576bac0..e929955 100755
--- a/plip/plipcmd
+++ b/plip/plipcmd
@@ -25,30 +25,29 @@ try:
     from plip.modules.preparation import *
     from plip.modules.visualize import visualize_in_pymol
     from plip.modules.plipremote import VisualizerData
-    from plip.modules.report import TextOutput
+    from plip.modules.report import StructureReport,__version__
     from plip.modules import config
     from plip.modules.mp import parallel_fn
+    from plip.modules.webservices import check_pdb_status, fetch_pdb
 except ImportError:
     from modules.preparation import *
-    from modules.visualize import visualize_in_pymol, PyMOLComplex
+    from modules.visualize import visualize_in_pymol
     from modules.plipremote import VisualizerData
-    from modules.report import TextOutput
+    from modules.report import StructureReport, __version__
     from modules import config
     from modules.mp import parallel_fn
+    from modules.webservices import check_pdb_status, fetch_pdb
 
 # Python standard library
 import sys
 import argparse
 from argparse import ArgumentParser
-import urllib2
 import time
 import multiprocessing
 import json
-
 # External libraries
 import lxml.etree as et
 
-__version__ = '1.3.1'
 descript = "Protein-Ligand Interaction Profiler (PLIP) v%s " \
            "is a command-line based tool to analyze interactions in a protein-ligand complex. " \
            "If you are using PLIP in your work, please cite: " \
@@ -63,41 +62,7 @@ def threshold_limiter(aparser, arg):
     return arg
 
 
-def check_pdb_status(pdbid):
-    """Returns the status and up-to-date entry in the PDB for a given PDB ID"""
-    url = 'http://www.rcsb.org/pdb/rest/idStatus?structureId=%s' % pdbid
-    xmlf = urllib2.urlopen(url)
-    xml = et.parse(xmlf)
-    xmlf.close()
-    status = None
-    current_pdbid = pdbid
-    for df in xml.xpath('//record'):
-        status = df.attrib['status']  # Status of an entry can be either 'UNKWOWN', 'OBSOLETE', or 'CURRENT'
-        if status == 'OBSOLETE':
-            current_pdbid = df.attrib['replacedBy']  # Contains the up-to-date PDB ID for obsolete entries
-    return [status, current_pdbid.lower()]
-
-
-def fetch_pdb(pdbid):
-    """Get the newest entry from the RCSB server for the given PDB ID. Exits with '1' if PDB ID is invalid."""
-    pdbid = pdbid.lower()
-    write_message('\nChecking status of PDB ID %s ... ' % pdbid)
-    state, current_entry = check_pdb_status(pdbid)  # Get state and current PDB ID
-
-    if state == 'OBSOLETE':
-        write_message('entry is obsolete, getting %s instead.\n' % current_entry)
-    elif state == 'CURRENT':
-        write_message('entry is up to date.\n')
-    elif state == 'UNKNOWN':
-        sysexit(3, 'Invalid PDB ID (Entry does not exist on PDB server)\n')
-    write_message('Downloading file from PDB ... ')
-    pdburl = 'http://www.rcsb.org/pdb/files/%s.pdb' % current_entry  # Get URL for current entry
-    pdbfile = urllib2.urlopen(pdburl).read()
-    # If no PDB file is available, a text is now shown with "We're sorry, but ..."
-    # Could previously be distinguished by an HTTP error
-    if 'sorry' in pdbfile:
-        sysexit(5, "No file in PDB format available from wwPDB for the given PDB ID.\n")
-    return [pdbfile, current_entry]
+
 
 
 def process_pdb(pdbfile, outpath):
@@ -111,50 +76,14 @@ def process_pdb(pdbfile, outpath):
     # #@todo Offers possibility for filter function from command line (by ligand chain, position, hetid)
     for ligand in mol.ligands:
         mol.characterize_complex(ligand)
-    excluded = mol.excluded
+
     create_folder_if_not_exists(outpath)
 
-    # Begin constructing the XML tree
-    report = et.Element('report')
-    plipversion = et.SubElement(report, 'plipversion')
-    plipversion.text = __version__
-    pdbid = et.SubElement(report, 'pdbid')
-    pdbid.text = mol.pymol_name.upper()
-    filetype = et.SubElement(report, 'filetype')
-    filetype.text = mol.filetype.upper()
-    pdbfile = et.SubElement(report, 'pdbfile')
-    pdbfile.text = mol.sourcefiles['pdbcomplex']
-    pdbfixes = et.SubElement(report, 'pdbfixes')
-    pdbfixes.text = str(mol.information['pdbfixes'])
-    filename = et.SubElement(report, 'filename')
-    filename.text = str(mol.sourcefiles['filename'])
-    exligs = et.SubElement(report, 'excluded_ligands')
-    for i, exlig in enumerate(excluded):
-        e = et.SubElement(exligs, 'excluded_ligand', id=str(i + 1))
-        e.text = exlig
-
-    # Write header of rST file
-    textlines = ['Prediction of noncovalent interactions for PDB structure %s' % mol.pymol_name.upper(), ]
-    textlines.append("=" * len(textlines[0]))
-    textlines.append('Created on %s using PLIP v%s\n' % (time.strftime("%Y/%m/%d"), __version__))
-    if len(excluded) != 0:
-        textlines.append('Excluded molecules as ligands: %s\n' % ','.join([lig for lig in excluded]))
+    # Generate the report files
+    streport = StructureReport(mol)
 
     config.MAXTHREADS = min(config.MAXTHREADS, len(mol.interaction_sets))
 
-    ####
-    # Create JSON dump of PyMOL helper class
-    ###
-
-    if config.VISJSON:
-        complexes = [VisualizerData(mol, site) for site in sorted(mol.interaction_sets)
-        if not len(mol.interaction_sets[site].interacting_res) == 0]
-
-        with open('%s/plipvis.json' % outpath, 'w') as f:
-            cplx_dict = {}
-            for cplx in complexes:
-                cplx_dict[cplx.uid] = cplx.to_json()
-            json.dump(cplx_dict, f)
 
     ######################################
     # PyMOL Visualization (parallelized) #
@@ -170,31 +99,11 @@ def process_pdb(pdbfile, outpath):
         else:
             [visualize_in_pymol(plcomplex) for plcomplex in complexes]
 
-    ##################################################################
-    # Generate XML- and rST-formatted reports for each binding site. #
-    ##################################################################
-
-    for i, site in enumerate(sorted(mol.interaction_sets)):
-        s = mol.interaction_sets[site]
-        bindingsite = TextOutput(s).generate_xml()
-        bindingsite.set('id', str(i + 1))
-        bindingsite.set('has_interactions', 'False')
-        report.insert(i + 1, bindingsite)
-        for itype in TextOutput(s).generate_txt():
-            textlines.append(itype)
-        if not s.no_interactions:
-            bindingsite.set('has_interactions', 'True')
-        else:
-            textlines.append('No interactions detected.')
-        sys.stdout = sys.__stdout__  # Change back to original stdout, gets changed when PyMOL has been used before
-
-    tree = et.ElementTree(report)
     if config.XML:  # Generate report in xml format
-        tree.write('%s/report.xml' % outpath, pretty_print=True, xml_declaration=True)
+        streport.write_xml()
 
     if config.TXT:  # Generate report in txt (rst) format
-        with open('%s/report.txt' % outpath, 'w') as f:
-            [f.write(textline + '\n') for textline in textlines]
+        streport.write_txt()
 
 def download_structure(inputpdbid):
     """Given a PDB ID, downloads the corresponding PDB structure.
@@ -297,6 +206,12 @@ if __name__ == '__main__':
     parser.add_argument("--nofix", dest="nofix", default=False,
                         help="Turns off fixing of PDB files.",
                         action="store_true")
+    parser.add_argument("--peptides", dest="peptides", default=[],
+                        help="Allows to define one or multiple chains as peptide ligands",
+                        nargs="+")
+    parser.add_argument("--keepmod", dest="keepmod", default=False,
+                        help="Keep modified residues as ligands",
+                        action="store_true")
     # Optional threshold arguments, not shown in help
     thr = namedtuple('threshold', 'name type')
     thresholds = [thr(name='aromatic_planarity', type='angle'),
@@ -327,7 +242,9 @@ if __name__ == '__main__':
     config.BASEPATH = config.OUTPATH  # Used for batch processing
     config.BREAKCOMPOSITE = arguments.breakcomposite
     config.ALTLOC = arguments.altlocation
+    config.PEPTIDES =arguments.peptides
     config.NOFIX = arguments.nofix
+    config.KEEPMOD = arguments.keepmod
     # Assign values to global thresholds
     for t in thresholds:
         tvalue = getattr(arguments, t.name)
diff --git a/plip/test/test_basic_functions.py b/plip/test/test_basic_functions.py
index 763cc47..79ffd51 100644
--- a/plip/test/test_basic_functions.py
+++ b/plip/test/test_basic_functions.py
@@ -17,10 +17,30 @@ See the License for the specific language governing permissions and
 limitations under the License.
 """
 
-
+# Python Standard Library
 import unittest
+import numpy
+import random
+
+# Own modules
 from plip.modules.preparation import PDBComplex
+from plip.modules.supplemental import euclidean3d, vector, vecangle, projection
+from plip.modules.supplemental import normalize_vector, cluster_doubles, centroid
 
+class TestLigandSupport(unittest.TestCase):
+    """Test for support of different ligands"""
+
+    def test_dna_rna(self):
+        """Test if DNA and RNA is correctly processed as ligands"""
+        tmpmol = PDBComplex()
+        tmpmol.load_pdb('./pdb/1tf6.pdb')
+        bsid = 'DA:B:1'
+        # DNA ligand four times consisting of 31 parts (composite)
+        self.assertEqual([len(ligand.members) for ligand in tmpmol.ligands].count(31), 4)
+        for ligset in [set((x[0] for x in ligand.members)) for ligand in tmpmol.ligands]:
+            if len(ligset) == 4:
+                # DNA only contains four bases
+                self.assertEqual(ligset, set(['DG', 'DC', 'DA', 'DT']))
 
 class TestMapping(unittest.TestCase):
     """Test"""
@@ -57,4 +77,70 @@ class TestMapping(unittest.TestCase):
                 self.assertEqual(contact.acc.o_orig_idx, 485)
             if contact.restype == 'LEU' and contact.resnr == 157:
                 self.assertEqual(contact.don.x_orig_idx, 1628)
-                self.assertEqual(contact.acc.o_orig_idx, 1191)
\ No newline at end of file
+                self.assertEqual(contact.acc.o_orig_idx, 1191)
+
+class GeometryTest(unittest.TestCase):
+    """Tests for geometrical calculations in PLIP"""
+
+    def vector_magnitude(self, v):
+        return numpy.sqrt(sum(x**2 for x in v))
+
+    # noinspection PyUnusedLocal
+    def setUp(self):
+        """Generate random data for the tests"""
+        # Generate two random n-dimensional float vectors, with -100 <= n <= 100 and values 0 <= i <= 1
+        dim = random.randint(1, 100)
+        self.rnd_vec = [random.uniform(-100, 100) for i in xrange(dim)]
+
+    def test_euclidean(self):
+        """Tests for mathematics.euclidean"""
+        # Are the results correct?
+        self.assertEqual(euclidean3d([0.0, 0.0, 0.0], [0.0, 0.0, 0.0]), 0)
+        self.assertEqual(euclidean3d([2.0, 3.0, 4.0], [2.0, 3.0, 4.0]), 0)
+        self.assertEqual(euclidean3d([4.0, 5.0, 6.0], [4.0, 5.0, 8.0]), 2.0)
+        # Does the function take vectors or tuples as an input? What about integers?
+        self.assertEqual(euclidean3d((4.0, 5.0, 6.0), [4.0, 5.0, 8.0]), 2.0)
+        self.assertEqual(euclidean3d((4.0, 5.0, 6.0), (4.0, 5.0, 8.0)), 2.0)
+        self.assertEqual(euclidean3d((4, 5, 6), (4.0, 5.0, 8.0)), 2.0)
+        # Is the output a float?
+        self.assertIsInstance(euclidean3d([2.0, 3.0, 4.0], [2.0, 3.0, 4.0]), float)
+
+    def test_vector(self):
+        """Tests for mathematics.vector"""
+        # Are the results correct?
+        self.assertEqual(list(vector([1, 1, 1], [0, 1, 0])), [-1, 0, -1])
+        self.assertEqual(list(vector([0, 0, 10], [0, 0, 4])), [0, 0, -6])
+        # Do I get an Numpy Array?
+        self.assertIsInstance(vector([1, 1, 1], [0, 1, 0]), numpy.ndarray)
+        # Do I get 'None' if the points have different dimensions?
+        self.assertEqual(vector([1, 1, 1], [0, 1, 0, 1]), None)
+
+    def test_vecangle(self):
+        """Tests for mathematics.vecangle"""
+        # Are the results correct?
+        self.assertEqual(vecangle([3, 4], [-8, 6], deg=False), numpy.radians(90.0))
+        self.assertEqual(vecangle([3, 4], [-8, 6]), 90.0)
+        self.assertAlmostEqual(vecangle([-1, -1], [1, 1], deg=False), numpy.pi)
+        # Correct if both vectors are equal?
+        self.assertEqual(vecangle([3, 3], [3, 3]), 0.0)
+
+    def test_centroid(self):
+        """Tests for mathematics.centroid"""
+        # Are the results correct?
+        self.assertEqual(centroid([[0, 0, 0], [2, 2, 2]]), [1.0, 1.0, 1.0])
+        self.assertEqual(centroid([[-5, 1, 2], [10, 2, 2]]), [2.5, 1.5, 2.0])
+
+    def test_normalize_vector(self):
+        """Tests for mathematics.normalize_vector"""
+        # Are the results correct?
+        self.assertAlmostEqual(self.vector_magnitude(normalize_vector(self.rnd_vec)), 1)
+
+    def test_projection(self):
+        """Tests for mathematics.projection"""
+        # Are the results correct?
+        self.assertEqual(projection([-1, 0, 0], [3, 3, 3], [1, 1, 1]), [3, 1, 1])
+
+    def test_cluster_doubles(self):
+        """Tests for mathematics.cluster_doubles"""
+        # Are the results correct?
+        self.assertEqual(set(cluster_doubles([(1, 3), (4, 1), (5, 6), (7, 5)])), {(1, 3, 4), (5, 6, 7)})
diff --git a/plip/test/test_literature_validated.py b/plip/test/test_literature_validated.py
index 1f10530..31c56a7 100644
--- a/plip/test/test_literature_validated.py
+++ b/plip/test/test_literature_validated.py
@@ -142,7 +142,7 @@ class LiteratureValidatedTest(unittest.TestCase):
         """
         tmpmol = PDBComplex()
         tmpmol.load_pdb('./pdb/2w0s.pdb')
-        bsid = 'MG:B:1206'  # Complex of BVDU with Magnesium Cofactor
+        bsid = 'BVP:B:1207'  # Complex of BVDU with Magnesium Cofactor
         for ligand in tmpmol.ligands:
             if ':'.join([ligand.hetid, ligand.chain, str(ligand.position)]) == bsid:
                 tmpmol.characterize_complex(ligand)
@@ -736,5 +736,3 @@ class LiteratureValidatedTest(unittest.TestCase):
         # Hydrogen bonds
         hbonds = {str(hbond.resnr)+hbond.reschain for hbond in s.hbonds_pdon+s.hbonds_ldon}
         self.assertTrue({'594A', '530A'}.issubset(hbonds))
-
-
diff --git a/plip/test/test_remote_services.py b/plip/test/test_remote_services.py
new file mode 100644
index 0000000..c8ecde8
--- /dev/null
+++ b/plip/test/test_remote_services.py
@@ -0,0 +1,46 @@
+# coding=utf-8
+"""
+Protein-Ligand Interaction Profiler - Analyze and visualize protein-ligand interactions in PDB files.
+test_remote_services.py - Unit Tests for remote services.
+Copyright 2014-2016 Sebastian Salentin
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+"""
+
+
+import unittest
+from plip.modules.webservices import check_pdb_status
+
+
+class TestPDB(unittest.TestCase):
+    """Test PDB Web Service methods"""
+
+    def test_pdb_entry_status(self):
+        # 1a0v is an obsolete entry and is replaced by 1y46
+        status, current_pdbid = check_pdb_status('1a0v')
+        self.assertEqual(status, 'OBSOLETE')
+        self.assertEqual(current_pdbid, '1y46')
+
+        # 1vsn is an current entry
+        status, current_pdbid = check_pdb_status('1vsn')
+        self.assertEqual(status, 'CURRENT')
+        self.assertEqual(current_pdbid, '1vsn')
+
+        # xxxx is not an PDB entry
+        status, current_pdbid = check_pdb_status('xxxx')
+        self.assertEqual(status, 'UNKNOWN')
+        self.assertEqual(current_pdbid, 'xxxx')
+
+
+        #print "Latest PDB ---", "passed" if get_latest_pdb("1a0v") == "1y46" else "not passed"
+#print "Resolution ---", "passed" if get_resolution("4HHB") == 1.74 else "not passed"
diff --git a/setup.py b/setup.py
index 1fff2dc..752d5cd 100644
--- a/setup.py
+++ b/setup.py
@@ -19,7 +19,7 @@ limitations under the License.
 from setuptools import setup
 
 setup(name='plip',
-      version='1.3.1',
+      version='1.3.2',
       description='PLIP - Fully automated protein-ligand interaction profiler',
       classifiers=[
           'Development Status :: 5 - Production/Stable',

-- 
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/plip.git



More information about the debian-med-commit mailing list