initial commit

5 files changed, 837 insertions, 0 deletions

diff --git a/plip/exchange/__init__.py b/plip/exchange/__init__.py
new file mode 100644
index 0000000..e69de29
--- /dev/null
+++ b/plip/exchange/__init__.py
diff --git a/plip/exchange/json.py b/plip/exchange/json.py
new file mode 100644
index 0000000..5c7ad5c
--- /dev/null
+++ b/plip/exchange/json.py
@@ -0,0 +1 @@
+# place holder for module to add Json support
\ No newline at end of file
diff --git a/plip/exchange/report.py b/plip/exchange/report.py
new file mode 100644
index 0000000..a559bf8
--- /dev/null
+++ b/plip/exchange/report.py
@@ -0,0 +1,497 @@
+import time
+from operator import itemgetter
+
+import lxml.etree as et
+
+from plip.basic import config
+from plip.basic.config import __version__
+from plip.structure.preparation import PDBComplex
+
+
+class StructureReport:
+    """Creates reports (xml or txt) for one structure/"""
+
+    def __init__(self, mol: PDBComplex, outputprefix: str = 'report'):
+        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
+        self.outputprefix = outputprefix
+
+    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"
+
+        maintainer_information = et.SubElement(report, 'maintainer_information')
+        maintainer_information.text = config.__maintainer__
+        mode = et.SubElement(report, 'mode')
+        if config.DNARECEPTOR:
+            mode.text = 'dna_receptor'
+        else:
+            mode.text = 'default'
+        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.get('filename') or None)
+        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
+        covalent = et.SubElement(report, 'covlinkages')
+        for i, covlinkage in enumerate(self.mol.covalent):
+            e = et.SubElement(covalent, 'covlinkage', id=str(i + 1))
+            f1 = et.SubElement(e, 'res1')
+            f2 = et.SubElement(e, 'res2')
+            f1.text = ":".join([covlinkage.id1, covlinkage.chain1, str(covlinkage.pos1)])
+            f2.text = ":".join([covlinkage.id2, covlinkage.chain2, str(covlinkage.pos2)])
+        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]))
+        if config.DNARECEPTOR:
+            textlines.append('DNA/RNA in structure was chosen as the receptor part.\n')
+        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.')
+
+    def write_xml(self, as_string=False):
+        """Write the XML report"""
+        if not as_string:
+            et.ElementTree(self.xmlreport).write('{}/{}.xml'.format(self.outpath, self.outputprefix), pretty_print=True,
+                                                 xml_declaration=True)
+        else:
+            output = et.tostring(self.xmlreport, pretty_print=True)
+            if config.RAWSTRING:
+                output = repr(output)
+            print(output)
+
+    def write_txt(self, as_string=False):
+        """Write the TXT report"""
+        if not as_string:
+            with open('{}/{}.txt'.format(self.outpath, self.outputprefix), 'w') as f:
+                [f.write(textline + '\n') for textline in self.txtreport]
+        else:
+            output = '\n'.join(self.txtreport)
+            if config.RAWSTRING:
+                output = repr(output)
+            print(output)
+
+
+class BindingSiteReport:
+    """Gather report data and generate reports for one binding site in different formats."""
+
+    def __init__(self, plcomplex):
+
+        ################
+        # GENERAL DATA #
+        ################
+
+        self.complex = plcomplex
+        self.ligand = self.complex.ligand
+        self.bindingsite = self.complex.bindingsite
+        self.output_path = self.complex.output_path
+        self.bsid = ':'.join([self.ligand.hetid, self.ligand.chain, str(self.ligand.position)])
+        self.longname = self.ligand.longname
+        self.ligtype = self.ligand.type
+        self.bs_res = self.bindingsite.bs_res
+        self.min_dist = self.bindingsite.min_dist
+        self.bs_res_interacting = self.complex.interacting_res
+        self.pdbid = self.complex.pdbid.upper()
+        self.lig_members = self.complex.lig_members
+        self.interacting_chains = self.complex.interacting_chains
+
+        ############################
+        # HYDROPHOBIC INTERACTIONS #
+        ############################
+
+        self.hydrophobic_features = (
+            'RESNR', 'RESTYPE', 'RESCHAIN', 'RESNR_LIG', 'RESTYPE_LIG', 'RESCHAIN_LIG', 'DIST', 'LIGCARBONIDX',
+            'PROTCARBONIDX', 'LIGCOO',
+            'PROTCOO')
+        self.hydrophobic_info = []
+        for hydroph in self.complex.hydrophobic_contacts:
+            self.hydrophobic_info.append((hydroph.resnr, hydroph.restype, hydroph.reschain, hydroph.resnr_l,
+                                          hydroph.restype_l, hydroph.reschain_l, '%.2f' % hydroph.distance,
+                                          hydroph.ligatom_orig_idx, hydroph.bsatom_orig_idx, hydroph.ligatom.coords,
+                                          hydroph.bsatom.coords))
+
+        ##################
+        # HYDROGEN BONDS #
+        ##################
+
+        self.hbond_features = (
+            'RESNR', 'RESTYPE', 'RESCHAIN', 'RESNR_LIG', 'RESTYPE_LIG', 'RESCHAIN_LIG', 'SIDECHAIN', 'DIST_H-A',
+            'DIST_D-A',
+            'DON_ANGLE',
+            'PROTISDON', 'DONORIDX', 'DONORTYPE', 'ACCEPTORIDX', 'ACCEPTORTYPE', 'LIGCOO', 'PROTCOO')
+        self.hbond_info = []
+        for hbond in self.complex.hbonds_pdon + self.complex.hbonds_ldon:
+            ligatom, protatom = (hbond.a, hbond.d) if hbond.protisdon else (hbond.d, hbond.a)
+            self.hbond_info.append((hbond.resnr, hbond.restype, hbond.reschain, hbond.resnr_l, hbond.restype_l,
+                                    hbond.reschain_l, hbond.sidechain,
+                                    '%.2f' % hbond.distance_ah, '%.2f' % hbond.distance_ad, '%.2f' % hbond.angle,
+                                    hbond.protisdon, hbond.d_orig_idx, hbond.dtype, hbond.a_orig_idx, hbond.atype,
+                                    ligatom.coords, protatom.coords))
+
+        #################
+        # WATER-BRIDGES #
+        #################
+
+        self.waterbridge_features = (
+            'RESNR', 'RESTYPE', 'RESCHAIN', 'RESNR_LIG', 'RESTYPE_LIG', 'RESCHAIN_LIG', 'DIST_A-W', 'DIST_D-W',
+            'DON_ANGLE',
+            'WATER_ANGLE',
+            'PROTISDON', 'DONOR_IDX', 'DONORTYPE', 'ACCEPTOR_IDX', 'ACCEPTORTYPE', 'WATER_IDX',
+            'LIGCOO', 'PROTCOO', 'WATERCOO')
+        # The coordinate format is an exception here, since the interaction is not only between ligand and protein
+        self.waterbridge_info = []
+        for wbridge in self.complex.water_bridges:
+            lig, prot = (wbridge.a, wbridge.d) if wbridge.protisdon else (wbridge.d, wbridge.a)
+            self.waterbridge_info.append((wbridge.resnr, wbridge.restype, wbridge.reschain, wbridge.resnr_l,
+                                          wbridge.restype_l, wbridge.reschain_l,
+                                          '%.2f' % wbridge.distance_aw, '%.2f' % wbridge.distance_dw,
+                                          '%.2f' % wbridge.d_angle, '%.2f' % wbridge.w_angle, wbridge.protisdon,
+                                          wbridge.d_orig_idx, wbridge.dtype, wbridge.a_orig_idx, wbridge.atype,
+                                          wbridge.water_orig_idx, lig.coords, prot.coords, wbridge.water.coords))
+
+        ################
+        # SALT BRIDGES #
+        ################
+
+        self.saltbridge_features = (
+            'RESNR', 'RESTYPE', 'RESCHAIN', 'RESNR_LIG', 'RESTYPE_LIG', 'RESCHAIN_LIG', 'DIST', 'PROTISPOS',
+            'LIG_GROUP',
+            'LIG_IDX_LIST',
+            'LIGCOO', 'PROTCOO')
+        self.saltbridge_info = []
+        for sb in self.complex.saltbridge_lneg + self.complex.saltbridge_pneg:
+            if sb.protispos:
+                group, ids = sb.negative.fgroup, [str(x) for x in sb.negative.atoms_orig_idx]
+                self.saltbridge_info.append((sb.resnr, sb.restype, sb.reschain, sb.resnr_l, sb.restype_l, sb.reschain_l,
+                                             '%.2f' % sb.distance, sb.protispos,
+                                             group.capitalize(), ",".join(ids),
+                                             tuple(sb.negative.center), tuple(sb.positive.center)))
+            else:
+                group, ids = sb.positive.fgroup, [str(x) for x in sb.positive.atoms_orig_idx]
+                self.saltbridge_info.append((sb.resnr, sb.restype, sb.reschain, sb.resnr_l, sb.restype_l, sb.reschain_l,
+                                             '%.2f' % sb.distance, sb.protispos,
+                                             group.capitalize(), ",".join(ids),
+                                             tuple(sb.positive.center), tuple(sb.negative.center)))
+
+        ###############
+        # PI-STACKING #
+        ###############
+
+        self.pistacking_features = (
+            'RESNR', 'RESTYPE', 'RESCHAIN', 'RESNR_LIG', 'RESTYPE_LIG', 'RESCHAIN_LIG', 'CENTDIST', 'ANGLE', 'OFFSET',
+            'TYPE',
+            'LIG_IDX_LIST', 'LIGCOO', 'PROTCOO')
+        self.pistacking_info = []
+        for stack in self.complex.pistacking:
+            ids = [str(x) for x in stack.ligandring.atoms_orig_idx]
+            self.pistacking_info.append((stack.resnr, stack.restype, stack.reschain, stack.resnr_l, stack.restype_l,
+                                         stack.reschain_l, '%.2f' % stack.distance,
+                                         '%.2f' % stack.angle, '%.2f' % stack.offset, stack.type, ",".join(ids),
+                                         tuple(stack.ligandring.center), tuple(stack.proteinring.center)))
+
+        ##########################
+        # PI-CATION INTERACTIONS #
+        ##########################
+
+        self.pication_features = (
+            'RESNR', 'RESTYPE', 'RESCHAIN', 'RESNR_LIG', 'RESTYPE_LIG', 'RESCHAIN_LIG', 'DIST', 'OFFSET', 'PROTCHARGED',
+            'LIG_GROUP',
+            'LIG_IDX_LIST', 'LIGCOO', 'PROTCOO')
+        self.pication_info = []
+        for picat in self.complex.pication_laro + self.complex.pication_paro:
+            if picat.protcharged:
+                ids = [str(x) for x in picat.ring.atoms_orig_idx]
+                group = 'Aromatic'
+                self.pication_info.append((picat.resnr, picat.restype, picat.reschain, picat.resnr_l, picat.restype_l,
+                                           picat.reschain_l, '%.2f' % picat.distance,
+                                           '%.2f' % picat.offset, picat.protcharged, group, ",".join(ids),
+                                           tuple(picat.ring.center), tuple(picat.charge.center)))
+            else:
+                ids = [str(x) for x in picat.charge.atoms_orig_idx]
+                group = picat.charge.fgroup
+                self.pication_info.append((picat.resnr, picat.restype, picat.reschain, picat.resnr_l, picat.restype_l,
+                                           picat.reschain_l, '%.2f' % picat.distance,
+                                           '%.2f' % picat.offset, picat.protcharged, group, ",".join(ids),
+                                           tuple(picat.charge.center), tuple(picat.ring.center)))
+
+        #################
+        # HALOGEN BONDS #
+        #################
+
+        self.halogen_features = (
+            'RESNR', 'RESTYPE', 'RESCHAIN', 'RESNR_LIG', 'RESTYPE_LIG', 'RESCHAIN_LIG', 'SIDECHAIN', 'DIST',
+            'DON_ANGLE',
+            'ACC_ANGLE',
+            'DON_IDX', 'DONORTYPE', 'ACC_IDX', 'ACCEPTORTYPE', 'LIGCOO', 'PROTCOO')
+        self.halogen_info = []
+        for halogen in self.complex.halogen_bonds:
+            self.halogen_info.append((halogen.resnr, halogen.restype, halogen.reschain, halogen.resnr_l,
+                                      halogen.restype_l, halogen.reschain_l, halogen.sidechain,
+                                      '%.2f' % halogen.distance, '%.2f' % halogen.don_angle, '%.2f' % halogen.acc_angle,
+                                      halogen.don_orig_idx, halogen.donortype,
+                                      halogen.acc_orig_idx, halogen.acctype,
+                                      halogen.acc.o.coords, halogen.don.x.coords))
+
+        ###################
+        # METAL COMPLEXES #
+        ###################
+
+        self.metal_features = (
+            'RESNR', 'RESTYPE', 'RESCHAIN', 'RESNR_LIG', 'RESTYPE_LIG', 'RESCHAIN_LIG', 'METAL_IDX', 'METAL_TYPE',
+            'TARGET_IDX', 'TARGET_TYPE',
+            'COORDINATION', 'DIST', 'LOCATION', 'RMS', 'GEOMETRY', 'COMPLEXNUM', 'METALCOO',
+            'TARGETCOO')
+        self.metal_info = []
+        # Coordinate format here is non-standard since the interaction partner can be either ligand or protein
+        for m in self.complex.metal_complexes:
+            self.metal_info.append(
+                (m.resnr, m.restype, m.reschain, m.resnr_l, m.restype_l, m.reschain_l, m.metal_orig_idx, m.metal_type,
+                 m.target_orig_idx, m.target_type, m.coordination_num, '%.2f' % m.distance,
+                 m.location, '%.2f' % m.rms, m.geometry, str(m.complexnum), m.metal.coords,
+                 m.target.atom.coords))
+
+    def write_section(self, name, features, info, f):
+        """Provides formatting for one section (e.g. hydrogen bonds)"""
+        if not len(info) == 0:
+            f.write('\n\n### %s ###\n' % name)
+            f.write('%s\n' % '\t'.join(features))
+            for line in info:
+                f.write('%s\n' % '\t'.join(map(str, line)))
+
+    def rst_table(self, array):
+        """Given an array, the function formats and returns and table in rST format."""
+        # Determine cell width for each column
+        cell_dict = {}
+        for i, row in enumerate(array):
+            for j, val in enumerate(row):
+                if j not in cell_dict:
+                    cell_dict[j] = []
+                cell_dict[j].append(val)
+        for item in cell_dict:
+            cell_dict[item] = max([len(x) for x in cell_dict[item]]) + 1  # Contains adapted width for each column
+
+        # Format top line
+        num_cols = len(array[0])
+        form = '+'
+        for col in range(num_cols):
+            form += (cell_dict[col] + 1) * '-'
+            form += '+'
+        form += '\n'
+
+        # Format values
+        for i, row in enumerate(array):
+            form += '| '
+            for j, val in enumerate(row):
+                cell_width = cell_dict[j]
+                form += str(val) + (cell_width - len(val)) * ' ' + '| '
+            form.rstrip()
+            form += '\n'
+
+            # Seperation lines
+            form += '+'
+            if i == 0:
+                sign = '='
+            else:
+                sign = '-'
+            for col in range(num_cols):
+                form += (cell_dict[col] + 1) * sign
+                form += '+'
+            form += '\n'
+        return form
+
+    def generate_txt(self):
+        """Generates an flat text report for a single binding site"""
+
+        txt = []
+        titletext = '%s (%s) - %s' % (self.bsid, self.longname, self.ligtype)
+        txt.append(titletext)
+        for i, member in enumerate(self.lig_members[1:]):
+            txt.append('  + %s' % ":".join(str(element) for element in member))
+        txt.append("-" * len(titletext))
+        txt.append("Interacting chain(s): %s\n" % ','.join([chain for chain in self.interacting_chains]))
+        for section in [['Hydrophobic Interactions', self.hydrophobic_features, self.hydrophobic_info],
+                        ['Hydrogen Bonds', self.hbond_features, self.hbond_info],
+                        ['Water Bridges', self.waterbridge_features, self.waterbridge_info],
+                        ['Salt Bridges', self.saltbridge_features, self.saltbridge_info],
+                        ['pi-Stacking', self.pistacking_features, self.pistacking_info],
+                        ['pi-Cation Interactions', self.pication_features, self.pication_info],
+                        ['Halogen Bonds', self.halogen_features, self.halogen_info],
+                        ['Metal Complexes', self.metal_features, self.metal_info]]:
+            iname, features, interaction_information = section
+            # Sort results first by res number, then by distance and finally ligand coordinates to get a unique order
+            interaction_information = sorted(interaction_information, key=itemgetter(0, 2, -2))
+            if not len(interaction_information) == 0:
+
+                txt.append('\n**%s**' % iname)
+                table = [features, ]
+                for single_contact in interaction_information:
+                    values = []
+                    for x in single_contact:
+                        if type(x) == str:
+                            values.append(x)
+                        elif type(x) == tuple and len(x) == 3:  # Coordinates
+                            values.append("%.3f, %.3f, %.3f" % x)
+                        else:
+                            values.append(str(x))
+                    table.append(values)
+                txt.append(self.rst_table(table))
+        txt.append('\n')
+        return txt
+
+    def generate_xml(self):
+        """Generates an XML-formatted report for a single binding site"""
+        report = et.Element('bindingsite')
+        identifiers = et.SubElement(report, 'identifiers')
+        longname = et.SubElement(identifiers, 'longname')
+        ligtype = et.SubElement(identifiers, 'ligtype')
+        hetid = et.SubElement(identifiers, 'hetid')
+        chain = et.SubElement(identifiers, 'chain')
+        position = et.SubElement(identifiers, 'position')
+        composite = et.SubElement(identifiers, 'composite')
+        members = et.SubElement(identifiers, 'members')
+        smiles = et.SubElement(identifiers, 'smiles')
+        inchikey = et.SubElement(identifiers, 'inchikey')
+
+        # Ligand properties. Number of (unpaired) functional atoms and rings.
+        lig_properties = et.SubElement(report, 'lig_properties')
+        num_heavy_atoms = et.SubElement(lig_properties, 'num_heavy_atoms')
+        num_hbd = et.SubElement(lig_properties, 'num_hbd')
+        num_hbd.text = str(self.ligand.num_hbd)
+        num_unpaired_hbd = et.SubElement(lig_properties, 'num_unpaired_hbd')
+        num_unpaired_hbd.text = str(self.complex.num_unpaired_hbd)
+        num_hba = et.SubElement(lig_properties, 'num_hba')
+        num_hba.text = str(self.ligand.num_hba)
+        num_unpaired_hba = et.SubElement(lig_properties, 'num_unpaired_hba')
+        num_unpaired_hba.text = str(self.complex.num_unpaired_hba)
+        num_hal = et.SubElement(lig_properties, 'num_hal')
+        num_hal.text = str(self.ligand.num_hal)
+        num_unpaired_hal = et.SubElement(lig_properties, 'num_unpaired_hal')
+        num_unpaired_hal.text = str(self.complex.num_unpaired_hal)
+        num_aromatic_rings = et.SubElement(lig_properties, 'num_aromatic_rings')
+        num_aromatic_rings.text = str(self.ligand.num_rings)
+        num_rot_bonds = et.SubElement(lig_properties, 'num_rotatable_bonds')
+        num_rot_bonds.text = str(self.ligand.num_rot_bonds)
+        molweight = et.SubElement(lig_properties, 'molweight')
+        molweight.text = str(self.ligand.molweight)
+        logp = et.SubElement(lig_properties, 'logp')
+        logp.text = str(self.ligand.logp)
+
+        ichains = et.SubElement(report, 'interacting_chains')
+        bsresidues = et.SubElement(report, 'bs_residues')
+        for i, ichain in enumerate(self.interacting_chains):
+            c = et.SubElement(ichains, 'interacting_chain', id=str(i + 1))
+            c.text = ichain
+        for i, bsres in enumerate(self.bs_res):
+            contact = 'True' if bsres in self.bs_res_interacting else 'False'
+            distance = '%.1f' % self.min_dist[bsres][0]
+            aatype = self.min_dist[bsres][1]
+            c = et.SubElement(bsresidues, 'bs_residue', id=str(i + 1), contact=contact, min_dist=distance, aa=aatype)
+            c.text = bsres
+        hetid.text, chain.text, position.text = self.ligand.hetid, self.ligand.chain, str(self.ligand.position)
+        composite.text = 'True' if len(self.lig_members) > 1 else 'False'
+        longname.text = self.longname
+        ligtype.text = self.ligtype
+        smiles.text = self.ligand.smiles
+        inchikey.text = self.ligand.inchikey
+        num_heavy_atoms.text = str(self.ligand.heavy_atoms)  # Number of heavy atoms in ligand
+        for i, member in enumerate(self.lig_members):
+            bsid = ":".join(str(element) for element in member)
+            m = et.SubElement(members, 'member', id=str(i + 1))
+            m.text = bsid
+        interactions = et.SubElement(report, 'interactions')
+
+        def format_interactions(element_name, features, interaction_information):
+            """Returns a formatted element with interaction information."""
+            interaction = et.Element(element_name)
+            # Sort results first by res number, then by distance and finally ligand coordinates to get a unique order
+            interaction_information = sorted(interaction_information, key=itemgetter(0, 2, -2))
+            for j, single_contact in enumerate(interaction_information):
+                if not element_name == 'metal_complexes':
+                    new_contact = et.SubElement(interaction, element_name[:-1], id=str(j + 1))
+                else:  # Metal Complex[es]
+                    new_contact = et.SubElement(interaction, element_name[:-2], id=str(j + 1))
+                for i, feature in enumerate(single_contact):
+                    # Just assign the value unless it's an atom list, use subelements in this case
+                    if features[i] == 'LIG_IDX_LIST':
+                        feat = et.SubElement(new_contact, features[i].lower())
+                        for k, atm_idx in enumerate(feature.split(',')):
+                            idx = et.SubElement(feat, 'idx', id=str(k + 1))
+                            idx.text = str(atm_idx)
+                    elif features[i].endswith('COO'):
+                        feat = et.SubElement(new_contact, features[i].lower())
+                        xc, yc, zc = feature
+                        xcoo = et.SubElement(feat, 'x')
+                        xcoo.text = '%.3f' % xc
+                        ycoo = et.SubElement(feat, 'y')
+                        ycoo.text = '%.3f' % yc
+                        zcoo = et.SubElement(feat, 'z')
+                        zcoo.text = '%.3f' % zc
+                    else:
+                        feat = et.SubElement(new_contact, features[i].lower())
+                        feat.text = str(feature)
+            return interaction
+
+        interactions.append(format_interactions('hydrophobic_interactions', self.hydrophobic_features,
+                                                self.hydrophobic_info))
+        interactions.append(format_interactions('hydrogen_bonds', self.hbond_features, self.hbond_info))
+        interactions.append(format_interactions('water_bridges', self.waterbridge_features, self.waterbridge_info))
+        interactions.append(format_interactions('salt_bridges', self.saltbridge_features, self.saltbridge_info))
+        interactions.append(format_interactions('pi_stacks', self.pistacking_features, self.pistacking_info))
+        interactions.append(format_interactions('pi_cation_interactions', self.pication_features, self.pication_info))
+        interactions.append(format_interactions('halogen_bonds', self.halogen_features, self.halogen_info))
+        interactions.append(format_interactions('metal_complexes', self.metal_features, self.metal_info))
+
+        # Mappings
+        mappings = et.SubElement(report, 'mappings')
+        smiles_to_pdb = et.SubElement(mappings, 'smiles_to_pdb')  # SMILES numbering to PDB file numbering (atoms)
+        bsid = ':'.join([self.ligand.hetid, self.ligand.chain, str(self.ligand.position)])
+        if self.ligand.atomorder is not None:
+            smiles_to_pdb_map = [(key, self.ligand.Mapper.mapid(self.ligand.can_to_pdb[key],
+                                                                mtype='protein', bsid=bsid)) for key in
+                                 self.ligand.can_to_pdb]
+            smiles_to_pdb.text = ','.join([str(mapping[0]) + ':' + str(mapping[1]) for mapping in smiles_to_pdb_map])
+        else:
+            smiles_to_pdb.text = ''
+
+        return report
diff --git a/plip/exchange/webservices.py b/plip/exchange/webservices.py
new file mode 100644
index 0000000..0c8cd3e
--- /dev/null
+++ b/plip/exchange/webservices.py
@@ -0,0 +1,54 @@
+import sys
+from urllib.error import HTTPError
+from urllib.request import urlopen
+
+import lxml.etree as et
+
+from plip.basic import logger
+
+logger = logger.get_logger()
+
+
+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 = 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()
+    logger.info(f'checking status of PDB-ID {pdbid}')
+    state, current_entry = check_pdb_status(pdbid)  # Get state and current PDB ID
+
+    if state == 'OBSOLETE':
+        logger.info(f'entry is obsolete, getting {current_entry} instead')
+    elif state == 'CURRENT':
+        logger.info('entry is up-to-date')
+    elif state == 'UNKNOWN':
+        logger.error('invalid PDB-ID (entry does not exist on PDB server)')
+        sys.exit(1)
+    logger.info('downloading file from PDB')
+    # get URL for current entry
+    # @todo needs update to react properly on response codes of RCSB servers
+    pdburl = f'http://www.rcsb.org/pdb/files/{current_entry}.pdb'
+    try:
+        pdbfile = urlopen(pdburl).read().decode()
+        # 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:
+            logger.error('no file in PDB format available from wwPDB for the given PDB ID.')
+            sys.exit(1)
+    except HTTPError:
+        logger.error('no file in PDB format available from wwPDB for the given PDB ID')
+        sys.exit(1)
+    return [pdbfile, current_entry]
diff --git a/plip/exchange/xml.py b/plip/exchange/xml.py
new file mode 100644
index 0000000..a5eee21
--- /dev/null
+++ b/plip/exchange/xml.py
@@ -0,0 +1,285 @@
+from lxml import etree
+
+from urllib.request import urlopen
+
+
+class XMLStorage:
+    """Generic class for storing XML data from PLIP XML files."""
+
+    def getdata(self, tree, location, force_string=False):
+        """Gets XML data from a specific element and handles types."""
+        found = tree.xpath('%s/text()' % location)
+        if not found:
+            return None
+        else:
+            data = found[0]
+        if force_string:
+            return data
+        if data == 'True':
+            return True
+        elif data == 'False':
+            return False
+        else:
+            try:
+                return int(data)
+            except ValueError:
+                try:
+                    return float(data)
+                except ValueError:
+                    # It's a string
+                    return data
+
+    def getcoordinates(self, tree, location):
+        """Gets coordinates from a specific element in PLIP XML"""
+        return tuple(float(x) for x in tree.xpath('.//%s/*/text()' % location))
+
+
+class Interaction(XMLStorage):
+    """Stores information on a specific interaction type"""
+
+    def __init__(self, interaction_part):
+        self.id = interaction_part.get('id')
+        self.resnr = self.getdata(interaction_part, 'resnr')
+        self.restype = self.getdata(interaction_part, 'restype', force_string=True)
+        self.reschain = self.getdata(interaction_part, 'reschain', force_string=True)
+        self.resnr_lig = self.getdata(interaction_part, 'resnr_lig')
+        self.restype_lig = self.getdata(interaction_part, 'restype_lig', force_string=True)
+        self.reschain_lig = self.getdata(interaction_part, 'reschain_lig', force_string=True)
+        self.ligcoo = self.getcoordinates(interaction_part, 'ligcoo')
+        self.protcoo = self.getcoordinates(interaction_part, 'protcoo')
+
+
+class HydrophobicInteraction(Interaction):
+    """Stores information on a hydrophobic interaction"""
+
+    def __init__(self, hydrophobic_part):
+        Interaction.__init__(self, hydrophobic_part)
+        self.dist = self.getdata(hydrophobic_part, 'dist')
+        self.ligcarbonidx = self.getdata(hydrophobic_part, 'ligcarbonidx')
+        self.protcarbonidx = self.getdata(hydrophobic_part, 'protcarbonidx')
+
+
+class HydrogenBond(Interaction):
+    """Stores information on a hydrogen bond interaction"""
+
+    def __init__(self, hbond_part):
+        Interaction.__init__(self, hbond_part)
+        self.sidechain = self.getdata(hbond_part, 'sidechain')
+        self.dist_h_a = self.getdata(hbond_part, 'dist_h-a')
+        self.dist_d_a = self.getdata(hbond_part, 'dist_d-a')
+        self.dist = self.dist_d_a
+
+        self.don_angle = self.getdata(hbond_part, 'don_angle')
+        self.protisdon = self.getdata(hbond_part, 'protisdon')
+        self.donoridx = self.getdata(hbond_part, 'donoridx')
+        self.acceptoridx = self.getdata(hbond_part, 'acceptoridx')
+        self.donortype = self.getdata(hbond_part, 'donortype', force_string=True)
+        self.acceptortype = self.getdata(hbond_part, 'acceptortype', force_string=True)
+
+
+class WaterBridge(Interaction):
+    """Stores information on a water bridge interaction"""
+
+    def __init__(self, wbridge_part):
+        Interaction.__init__(self, wbridge_part)
+        self.dist_a_w = self.getdata(wbridge_part, 'dist_a-w')
+        self.dist_d_w = self.getdata(wbridge_part, 'dist_d-w')
+        self.don_angle = self.getdata(wbridge_part, 'don_angle')
+        self.water_angle = self.getdata(wbridge_part, 'water_angle')
+        self.protisdon = self.getdata(wbridge_part, 'protisdon')
+        self.dist = self.dist_a_w if self.protisdon else self.dist_d_w
+
+        self.donor_idx = self.getdata(wbridge_part, 'donor_idx')
+        self.acceptor_idx = self.getdata(wbridge_part, 'acceptor_idx')
+        self.donortype = self.getdata(wbridge_part, 'donortype', force_string=True)
+        self.acceptortype = self.getdata(wbridge_part, 'acceptortype', force_string=True)
+        self.water_idx = self.getdata(wbridge_part, 'water_idx')
+        self.watercoo = self.getcoordinates(wbridge_part, 'watercoo')
+
+
+class SaltBridge(Interaction):
+    """Stores information on a salt bridge interaction"""
+
+    def __init__(self, sbridge_part):
+        Interaction.__init__(self, sbridge_part)
+        self.dist = self.getdata(sbridge_part, 'dist')
+        self.protispos = self.getdata(sbridge_part, 'protispos')
+        self.lig_group = self.getdata(sbridge_part, 'lig_group', force_string=True)
+        self.lig_idx_list = [int(tagpart.text) for tagpart in
+                             sbridge_part.xpath('lig_idx_list/idx')]
+
+
+class PiStacking(Interaction):
+    """Stores information on a pi stacking interaction"""
+
+    def __init__(self, pistack_part):
+        Interaction.__init__(self, pistack_part)
+        self.centdist = self.getdata(pistack_part, 'centdist')
+        self.dist = self.centdist
+        self.angle = self.getdata(pistack_part, 'angle')
+        self.offset = self.getdata(pistack_part, 'offset')
+        self.type = self.getdata(pistack_part, 'type')
+        self.lig_idx_list = [int(tagpart.text) for tagpart in
+                             pistack_part.xpath('lig_idx_list/idx')]
+
+
+class PiCation(Interaction):
+    """Stores information on a pi cation interaction"""
+
+    def __init__(self, pication_part):
+        Interaction.__init__(self, pication_part)
+        self.dist = self.getdata(pication_part, 'dist')
+        self.offset = self.getdata(pication_part, 'offset')
+        self.protcharged = self.getdata(pication_part, 'protcharged')
+        self.lig_group = self.getdata(pication_part, 'lig_group')
+        self.lig_idx_list = [int(tag.text) for tag in pication_part.xpath('.//lig_idx_list/idx')]
+
+
+class HalogenBond(Interaction):
+    """Stores information on a halogen bond interaction"""
+
+    def __init__(self, halogen_part):
+        Interaction.__init__(self, halogen_part)
+        self.dist = self.getdata(halogen_part, 'dist')
+        self.don_angle = self.getdata(halogen_part, 'don_angle')
+        self.acc_angle = self.getdata(halogen_part, 'acc_angle')
+        self.donortype = self.getdata(halogen_part, 'donortype', force_string=True)
+        self.acceptortype = self.getdata(halogen_part, 'acceptortype', force_string=True)
+        self.don_idx = self.getdata(halogen_part, 'don_idx')
+        self.acc_idx = self.getdata(halogen_part, 'acc_idx')
+        self.sidechain = self.getdata(halogen_part, 'sidechain')
+
+
+class MetalComplex(Interaction):
+    """Stores information on a metal complexe interaction"""
+
+    def __init__(self, metalcomplex_part):
+        Interaction.__init__(self, metalcomplex_part)
+        self.metal_idx = self.getdata(metalcomplex_part, 'metal_idx')
+        self.metal_type = self.getdata(metalcomplex_part, 'metal_type', force_string=True)
+        self.target_idx = self.getdata(metalcomplex_part, 'target_idx')
+        self.target_type = self.getdata(metalcomplex_part, 'target_type', force_string=True)
+        self.coordination = self.getdata(metalcomplex_part, 'coordination')
+        self.dist = self.getdata(metalcomplex_part, 'dist')
+        self.location = self.getdata(metalcomplex_part, 'location', force_string=True)
+        self.rms = self.getdata(metalcomplex_part, 'rms')
+        self.geometry = self.getdata(metalcomplex_part, 'geometry', force_string=True)
+        self.complexnum = self.getdata(metalcomplex_part, 'complexnum')
+        self.targetcoo = self.getcoordinates(metalcomplex_part, 'targetcoo')
+        self.metalcoo = self.getcoordinates(metalcomplex_part, 'metalcoo')
+
+
+class BSite(XMLStorage):
+    """Stores all information about an specific binding site."""
+
+    def __init__(self, bindingsite, pdbid):
+        self.bindingsite = bindingsite
+        self.pdbid = pdbid
+        self.bsid = ":".join(bindingsite.xpath('identifiers/*/text()')[2:5])
+        self.uniqueid = ":".join([self.pdbid, self.bsid])
+        self.hetid = self.getdata(bindingsite, 'identifiers/hetid', force_string=True)
+        self.longname = self.getdata(bindingsite, 'identifiers/longname', force_string=True)
+        self.ligtype = self.getdata(bindingsite, 'identifiers/ligtype', force_string=True)
+        self.smiles = self.getdata(bindingsite, 'identifiers/smiles', force_string=True)
+        self.inchikey = self.getdata(bindingsite, 'identifiers/inchikey', force_string=True)
+        self.position = self.getdata(bindingsite, 'identifiers/position')
+        self.chain = self.getdata(bindingsite, 'identifiers/chain', force_string=True)
+
+        # Information on binding site members
+        self.members = []
+        for member in bindingsite.xpath('identifiers/members/member'):
+            self.members += member.xpath('text()')
+
+        self.composite = self.getdata(bindingsite, 'identifiers/composite')
+
+        # Ligand Properties
+        self.heavy_atoms = self.getdata(bindingsite, 'lig_properties/num_heavy_atoms')
+        self.hbd = self.getdata(bindingsite, 'lig_properties/num_hbd')
+        self.unpaired_hbd = self.getdata(bindingsite, 'lig_properties/num_unpaired_hbd')
+        self.hba = self.getdata(bindingsite, 'lig_properties/num_hba')
+        self.unpaired_hba = self.getdata(bindingsite, 'lig_properties/num_unpaired_hba')
+        self.hal = self.getdata(bindingsite, 'lig_properties/num_hal')
+        self.unpaired_hal = self.getdata(bindingsite, 'lig_properties/num_unpaired_hal')
+        self.molweight = self.getdata(bindingsite, 'lig_properties/molweight')
+        self.logp = self.getdata(bindingsite, 'lig_properties/logp')
+        self.rotatable_bonds = self.getdata(bindingsite, 'lig_properties/num_rotatable_bonds')
+        self.rings = self.getdata(bindingsite, 'lig_properties/num_aromatic_rings')
+
+        # Binding Site residues
+        self.bs_res = []
+        for tagpart in bindingsite.xpath('bs_residues/bs_residue'):
+            resnumber, reschain = tagpart.text[:-1], tagpart.text[-1]
+            aa, contact, min_dist = tagpart.get('aa'), tagpart.get('contact'), tagpart.get('min_dist')
+            new_bs_res = {'resnr': int(resnumber), 'reschain': reschain, 'aa': aa,
+                          'contact': True if contact == 'True' else False, 'min_dist': float(min_dist)}
+            self.bs_res.append(new_bs_res)
+
+        # Interacting chains
+        self.interacting_chains = []
+        for chain in bindingsite.xpath('interacting_chains/interacting_chain'):
+            self.interacting_chains += chain.xpath('text()')
+
+        # Interactions
+        interactions = bindingsite.xpath('interactions')[0]
+        self.hydrophobics = [HydrophobicInteraction(x) for x in
+                             interactions.xpath('hydrophobic_interactions/hydrophobic_interaction')]
+        self.hbonds = [HydrogenBond(x) for x in interactions.xpath('hydrogen_bonds/hydrogen_bond')]
+        self.wbridges = [WaterBridge(x) for x in interactions.xpath('water_bridges/water_bridge')]
+        self.sbridges = [SaltBridge(x) for x in interactions.xpath('salt_bridges/salt_bridge')]
+        self.pi_stacks = [PiStacking(x) for x in interactions.xpath('pi_stacks/pi_stack')]
+        self.pi_cations = [PiCation(x) for x in interactions.xpath('pi_cation_interactions/pi_cation_interaction')]
+        self.halogens = [HalogenBond(x) for x in interactions.xpath('halogen_bonds/halogen_bond')]
+        self.metal_complexes = [MetalComplex(x) for x in interactions.xpath('metal_complexes/metal_complex')]
+        self.num_contacts = len(self.hydrophobics) + len(self.hbonds) + len(self.wbridges) + len(self.sbridges) + \
+            len(self.pi_stacks) + len(self.pi_cations) + len(self.halogens) + len(self.metal_complexes)
+        self.has_interactions = self.num_contacts > 0
+
+        self.get_atom_mapping()
+        self.counts = self.get_counts()
+
+    def get_atom_mapping(self):
+        """Parses the ligand atom mapping."""
+        # Atom mappings
+        smiles_to_pdb_mapping = self.bindingsite.xpath('mappings/smiles_to_pdb/text()')
+        if not smiles_to_pdb_mapping:
+            self.mappings = {'smiles_to_pdb': None, 'pdb_to_smiles': None}
+        else:
+            smiles_to_pdb_mapping = {int(y[0]): int(y[1]) for y in [x.split(':')
+                                                                    for x in smiles_to_pdb_mapping[0].split(',')]}
+            self.mappings = {'smiles_to_pdb': smiles_to_pdb_mapping}
+            self.mappings['pdb_to_smiles'] = {v: k for k, v in self.mappings['smiles_to_pdb'].items()}
+
+    def get_counts(self):
+        """counts the interaction types and backbone hydrogen bonding in a binding site"""
+
+        hbondsback = len([hb for hb in self.hbonds if not hb.sidechain])
+        counts = {'hydrophobics': len(self.hydrophobics), 'hbonds': len(self.hbonds),
+                  'wbridges': len(self.wbridges), 'sbridges': len(self.sbridges), 'pistacks': len(self.pi_stacks),
+                  'pications': len(self.pi_cations), 'halogens': len(self.halogens), 'metal': len(self.metal_complexes),
+                  'hbond_back': hbondsback, 'hbond_nonback': (len(self.hbonds) - hbondsback)}
+        counts['total'] = counts['hydrophobics'] + counts['hbonds'] + counts['wbridges'] + \
+            counts['sbridges'] + counts['pistacks'] + counts['pications'] + counts['halogens'] + counts['metal']
+        return counts
+
+
+class PlipXML(XMLStorage):
+    """Parses and stores all information from a PLIP XML file."""
+
+    def __init__(self, xmlfile):
+        self.load_data(xmlfile)
+
+        # Parse general information
+        self.version = self.getdata(self.doc, '/report/plipversion/')
+        self.pdbid = self.getdata(self.doc, '/report/pdbid', force_string=True)
+        self.filetype = self.getdata(self.doc, '/report/filetype')
+        self.fixed = self.getdata(self.doc, '/report/pdbfixes/')
+        self.filename = self.getdata(self.doc, '/report/filename')
+        self.excluded = self.doc.xpath('/report/excluded_ligands/excluded_ligand/text()')
+
+        # Parse binding site information
+        self.bsites = {BSite(bs, self.pdbid).bsid: BSite(bs, self.pdbid) for bs in self.doc.xpath('//bindingsite')}
+        self.num_bsites = len(self.bsites)
+
+    def load_data(self, xmlfile):
+        """Loads/parses an XML file and saves it as a tree if successful."""
+        self.doc = etree.parse(xmlfile)
\ No newline at end of file