view @ 14:c19608db6b15 draft default tip

"planemo upload for repository commit 327c29cc43f56d7067ab9fa51323ea31951db98b"
author bgruening
date Tue, 10 Nov 2020 20:38:46 +0000
parents bd678d7db2ae
line wrap: on
line source

#!/usr/bin/env python
    Small library with cheminformatic functions based on openbabel and pgchem.
    Copyright 2012, Bjoern Gruening and Xavier Lucas

import glob
import re
import subprocess
import sys
import tempfile
from multiprocessing import Pool

    from galaxy import eggs
except ImportError:
    psycopg2 = None
    print('psycopg2 is not available. It is currently used in the pgchem wrappers, that are not shipped with default CTB')

    from openbabel import openbabel, pybel
except ImportError:
    openbabel, pybel = None, None
    print('OpenBabel could not be found. A few functions are not available without OpenBabel.')

def CountLines(path):
    out = subprocess.Popen(['wc', '-l', path],
    return int(out.partition(b' ')[0])

def grep(pattern, file_obj):
    grepper = re.compile(pattern)
    for line in file_obj:
            return True
    return False

def check_filetype(filepath):
    mol = False
    possible_inchi = True
    for line_counter, line in enumerate(open(filepath)):
        if line_counter > 10000:
        if line.find('$$$$') != -1:
            return 'sdf'
        elif line.find('@<TRIPOS>MOLECULE') != -1:
            return 'mol2'
        elif line.find('ligand id') != -1:
            return 'drf'
        elif possible_inchi and re.findall('^InChI=', line):
            return 'inchi'
        elif re.findall(r'^M\s+END', line):
            mol = True
        # first line is not an InChI, so it can't be an InChI file
        possible_inchi = False

    if mol:
        # END can occures before $$$$, so and SDF file will
        # be recognised as mol, if you not using this hack'
        return 'mol'
    return 'smi'

def db_connect(args):
        db_conn = psycopg2.connect("dbname=%s user=%s host=%s password=%s" % (args.dbname, args.dbuser, args.dbhost, args.dbpasswd))
        return db_conn
    except psycopg2.Error:
        sys.exit('Unable to connect to the db')

ColumnNames = {
    'can_smiles': 'Canonical SMILES',
    'can': 'Canonical SMILES',
    'inchi': 'InChI',
    'inchi_key': 'InChI key',
    'inchi_key_first': 'InChI key first',
    'inchi_key_last': 'InChI key last',
    'molwt': 'Molecular weight',
    'hbd': 'Hydrogen-bond donors',
    'donors': 'Hydrogen-bond donors',
    'hba': 'Hydrogen-bond acceptors',
    'acceptors': 'Hydrogen-bond acceptors',
    'rotbonds': 'Rotatable bonds',
    'logp': 'logP',
    'psa': 'Polar surface area',
    'mr': 'Molecular refractivity',
    'atoms': 'Number of heavy atoms',
    'rings': 'Number of rings',
    'set_bits': 'FP2 bits',
    'id': 'Internal identifier',
    'tani': 'Tanimoto coefficient',
    'spectrophore': 'Spectrophores(TM)',
    'dist_spectrophore': 'Spectrophores(TM) distance to target',
    'synonym': 'Entry id',

OBDescriptor = {
    'atoms': ["atoms", "Number of atoms"],
    'hatoms': ["hatoms", "Number of heavy atoms"],  # self defined tag hatoms in plugindefines.txt
    'can_smiles': ["cansmi", "Canonical SMILES"],
    'can_smilesNS': ["cansmiNS", "Canonical SMILES without isotopes or stereo"],
    # ["abonds", "Number of aromatic bonds"],
    # ["bonds", "Number of bonds"],
    # ["dbonds", "Number of double bonds"],
    # ["formula", "Chemical formula"],
    'hba': ["HBA1", "Number of Hydrogen Bond Acceptors 1 (JoelLib)"],
    'hba2': ["HBA2", "Number of Hydrogen Bond Acceptors 2 (JoelLib)"],
    'hbd': ["HBD", "Number of Hydrogen Bond Donors (JoelLib)"],
    'inchi': ["InChI", "IUPAC InChI identifier"],
    'inchi_key': ["InChIKey", "InChIKey"],
    # ["L5", "Lipinski Rule of Five"],
    'logp': ["logP", "octanol/water partition coefficient"],
    'mr': ["MR", "molar refractivity"],
    'molwt': ["MW", "Molecular Weight filter"],
    # ["nF", "Number of Fluorine Atoms"],
    # ["s", "SMARTS filter"],
    # ["sbonds", "Number of single bonds"],
    # ["smarts", "SMARTS filter"],
    # ["tbonds", "Number of triple bonds"],
    # ["title", "For comparing a molecule's title"],
    'psa': ["TPSA", "topological polar surface area"],
    'rotbonds': ['ROTATABLE_BOND', 'rotatable bonds'],

def print_output(args, rows):
    if args.oformat == 'table':
        outfile = open(args.output, 'w')
        requested_fields = (filter(lambda x: x not in ["[", "]", "'"], args.fetch)).split(', ')
        if args.header:
            outfile.write('Identifier\t' + '\t'.join([ColumnNames[key] for key in requested_fields]) + '\n')
        for row in rows:
            outfile.write(row['synonym'] + '\t' + '\t'.join([str(row[key]) for key in requested_fields]) + '\n')

    elif args.oformat in ['sdf', 'mol2']:
        outfile = pybel.Outputfile(args.oformat, args.output, overwrite=True)
        for row in rows:
                mol = pybel.readstring('sdf', row['mol'])
                if args.oformat == 'sdf':
                    keys = filter(lambda x: x not in ["[", "]", "'"], args.fetch).split(', ')
          {ColumnNames['synonym']: row['synonym']})
                    if 'inchi_key' in keys:
                        keys = (', '.join(keys).replace("inchi_key", "inchi_key_first, inchi_key_last")).split(', ')
                    [{ColumnNames[key]: row[key]}) for key in keys if key]
            except OSError:
        outfile = open(args.output, 'w')
        outfile.write('\n'.join(['%s\t%s' % (row[args.oformat], row['synonym']) for row in rows]))

def pybel_stop_logging():

def get_properties_ext(mol):
    HBD = pybel.Smarts("[!#6;!H0]")
    HBA = pybel.Smarts(("[$([$([#8,#16]);!$(*=N~O);"
    calc_desc_dict = mol.calcdesc()

        logp = calc_desc_dict['logP']
    except KeyError:
        logp = calc_desc_dict['LogP']

    return {"molwt": mol.molwt,
            "logp": logp,
            "donors": len(HBD.findall(mol)),
            "acceptors": len(HBA.findall(mol)),
            "psa": calc_desc_dict['TPSA'],
            "mr": calc_desc_dict['MR'],
            "rotbonds": mol.OBMol.NumRotors(),
            "can": mol.write("can").split()[0].strip(),  # tthis one works fine for both zinc and chembl (no ZINC code added after can descriptor string)
            "inchi": mol.write("inchi").strip(),
            "inchi_key": get_inchikey(mol).strip(),
            "rings": len(mol.sssr),
            "atoms": mol.OBMol.NumHvyAtoms(),
            "spectrophore": OBspectrophore(mol),

def get_inchikey(mol):
    conv = openbabel.OBConversion()
    conv.SetInAndOutFormats("mol", "inchi")
    conv.SetOptions("K", conv.OUTOPTIONS)
    inchikey = conv.WriteString(mol.OBMol)
    return inchikey

def OBspectrophore(mol):
    spectrophore = pybel.ob.OBSpectrophore()
    # Parameters: rotation angle = 20, normalization for mean and sd, accuracy = 3.0 A and non-stereospecific cages.
    return ', '.join(["%.3f" % value for value in spectrophore.GetSpectrophore(mol.OBMol)])

def split_library(lib_path, lib_format='sdf', package_size=None):
    Split a library of compounds. Usage: split_library(lib_path, lib_format, package_size)
    IT currently ONLY WORKS FOR SD-Files
    pack = 1
    mol_counter = 0

    outfile = open('/%s/%s_pack_%i.%s' % ('/'.join(lib_path.split('/')[:-1]), lib_path.split('/')[-1].split('.')[0], pack, 'sdf'), 'w')

    for line in open(lib_path, 'r'):
        if line.strip() == '$$$$':
            mol_counter += 1
            if mol_counter % package_size == 0:
                pack += 1
                outfile = open('/%s/%s_pack_%i.%s' % ('/'.join(lib_path.split('/')[:-1]), lib_path.split('/')[-1].split('.')[0], pack, 'sdf'), 'w')
                if mol_counter * 10 % package_size == 0:
                    print('%i molecules parsed, starting pack nr. %i' % (mol_counter, pack - 1))

    return True

def split_smi_library(smiles_file, structures_in_one_file):
    Split a file with SMILES to several files for multiprocessing usage.
    Usage: split_smi_library(smiles_file, 10)
    output_files = []
    tfile = tempfile.NamedTemporaryFile(delete=False)

    smiles_handle = open(smiles_file, 'r')
    for count, line in enumerate(smiles_handle):
        if count % structures_in_one_file == 0 and count != 0:
            tfile = tempfile.NamedTemporaryFile(delete=False)
    return output_files

def mp_run(input_path, regex, PROCESSES, function_to_call):
    paths = []
    [paths.append(compound_file) for compound_file in glob.glob(str(input_path) + str(regex))]

    pool = Pool(processes=PROCESSES)
    print('Process initialized with', PROCESSES, 'processors')
    result = pool.map_async(function_to_call, paths)

    return paths

if __name__ == '__main__':