cheminfolib.py

#!/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

import numpy as np
try:
    from galaxy import eggs

    eggs.require("psycopg2")
except Exception:
    psycopg2 = None
    print(
        "psycopg2 is not available. It is currently used in the pgchem wrappers, that are not shipped with default CTB"
    )

try:
    from openbabel import openbabel, pybel
except Exception:
    print(
        "OpenBabel could not be found. A few functions are not available without OpenBabel."
    )


def CountLines(path):
    out = subprocess.Popen(
        ["wc", "-l", path], stdout=subprocess.PIPE, stderr=subprocess.STDOUT
    ).communicate()[0]
    return int(out.partition(b" ")[0])


def grep(pattern, file_obj):
    grepper = re.compile(pattern)
    for line in file_obj:
        if grepper.search(line):
            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:
            break
        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("^M\s+END", line):  # noqa W605
            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):
    try:
        db_conn = psycopg2.connect(
            "dbname=%s user=%s host=%s password=%s"
            % (args.dbname, args.dbuser, args.dbhost, args.dbpasswd)
        )
        return db_conn
    except Exception:
        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:
            try:
                mol = pybel.readstring("sdf", row["mol"])
                if args.oformat == "sdf":
                    keys = filter(lambda x: x not in ["[", "]", "'"], args.fetch).split(
                        ", "
                    )
                    mol.data.update({ColumnNames["synonym"]: row["synonym"]})
                    if "inchi_key" in keys:
                        keys = (
                            ", ".join(keys).replace(
                                "inchi_key", "inchi_key_first, inchi_key_last"
                            )
                        ).split(", ")
                    [
                        mol.data.update({ColumnNames[key]: row[key]})
                        for key in keys
                        if key
                    ]
                outfile.write(mol)
            except Exception:
                pass
    else:
        outfile = open(args.output, "w")
        outfile.write(
            "\n".join(["%s\t%s" % (row[args.oformat], row["synonym"]) for row in rows])
        )
    outfile.close()


def pybel_stop_logging():
    openbabel.obErrorLog.StopLogging()


def get_properties_ext(mol):

    HBD = pybel.Smarts("[!#6;!H0]")
    HBA = pybel.Smarts(
        "[$([$([#8,#16]);!$(*=N~O);"
        + "!$(*~N=O);X1,X2]),$([#7;v3;"
        + "!$([nH]);!$(*(-a)-a)])]"
    )
    calc_desc_dict = mol.calcdesc()

    try:
        logp = calc_desc_dict["logP"]
    except Exception:
        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(),  # this 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.
    spectrophore.SetNormalization(spectrophore.NormalizationTowardsZeroMeanAndUnitStd)
    return ", ".join(
        ["%.3f" % value for value in spectrophore.GetSpectrophore(mol.OBMol)]
    )


def squared_euclidean_distance(a, b):
    try:
        return ((np.asarray(a) - np.asarray(b)) ** 2).sum()
    except ValueError:
        return 0


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"):
        outfile.write(line)
        if line.strip() == "$$$$":
            mol_counter += 1
            if mol_counter % package_size == 0:
                outfile.close()
                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)
                    )
    outfile.close()

    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.close()
            output_files.append(tfile.name)
            tfile = tempfile.NamedTemporaryFile(delete=False)
        tfile.write(bytes(line, "utf-8"))
    tfile.close()
    output_files.append(tfile.name)
    smiles_handle.close()
    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))
    ]
    paths.sort()

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

    return paths


if __name__ == "__main__":
    print(check_filetype(sys.argv[1]))