annotate qed.py @ 3:52a8d34dd08f draft default tip

"planemo upload for repository https://github.com/bgruening/galaxytools/tree/master/chemicaltoolbox/silicos-it/qed commit 943ff93be2257426d69a8406ed55c838495ecf3f"

#!/usr/bin/env python
__all__ = ["weights_max", "weights_mean", "weights_none", "default"]

import argparse
import os
import re
import sys
# General
from copy import deepcopy
from math import exp, log

from rdkit import Chem
# RDKit
from rdkit.Chem import Descriptors


class SilicosItError(Exception):
    """Base class for exceptions in Silicos-it code"""


class WrongArgument(SilicosItError):
    """
    Exception raised when argument to function is not of correct type.

    Attributes:
        function -- function in which error occurred
        msg      -- explanation of the error
    """

    def __init__(self, function, msg):
        self.function = function
        self.msg = msg


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"


AliphaticRings = Chem.MolFromSmarts("[$([A;R][!a])]")

AcceptorSmarts = [
    "[oH0;X2]",
    "[OH1;X2;v2]",
    "[OH0;X2;v2]",
    "[OH0;X1;v2]",
    "[O-;X1]",
    "[SH0;X2;v2]",
    "[SH0;X1;v2]",
    "[S-;X1]",
    "[nH0;X2]",
    "[NH0;X1;v3]",
    "[$([N;+0;X3;v3]);!$(N[C,S]=O)]",
]
Acceptors = []
for hba in AcceptorSmarts:
    Acceptors.append(Chem.MolFromSmarts(hba))

StructuralAlertSmarts = [
    "*1[O,S,N]*1",
    "[S,C](=[O,S])[F,Br,Cl,I]",
    "[CX4][Cl,Br,I]",
    "[C,c]S(=O)(=O)O[C,c]",
    "[$([CH]),$(CC)]#CC(=O)[C,c]",
    "[$([CH]),$(CC)]#CC(=O)O[C,c]",
    "n[OH]",
    "[$([CH]),$(CC)]#CS(=O)(=O)[C,c]",
    "C=C(C=O)C=O",
    "n1c([F,Cl,Br,I])cccc1",
    "[CH1](=O)",
    "[O,o][O,o]",
    "[C;!R]=[N;!R]",
    "[N!R]=[N!R]",
    "[#6](=O)[#6](=O)",
    "[S,s][S,s]",
    "[N,n][NH2]",
    "C(=O)N[NH2]",
    "[C,c]=S",
    "[$([CH2]),$([CH][CX4]),$(C([CX4])[CX4])]=[$([CH2]),$([CH][CX4]),$(C([CX4])[CX4])]",
    "C1(=[O,N])C=CC(=[O,N])C=C1",
    "C1(=[O,N])C(=[O,N])C=CC=C1",
    "a21aa3a(aa1aaaa2)aaaa3",
    "a31a(a2a(aa1)aaaa2)aaaa3",
    "a1aa2a3a(a1)A=AA=A3=AA=A2",
    "c1cc([NH2])ccc1",
    "[Hg,Fe,As,Sb,Zn,Se,se,Te,B,Si,Na,Ca,Ge,Ag,Mg,K,Ba,Sr,Be,Ti,Mo,Mn,Ru,Pd,Ni,Cu,Au,Cd,Al,Ga,Sn,Rh,Tl,Bi,Nb,Li,Pb,Hf,Ho]",
    "I",
    "OS(=O)(=O)[O-]",
    "[N+](=O)[O-]",
    "C(=O)N[OH]",
    "C1NC(=O)NC(=O)1",
    "[SH]",
    "[S-]",
    "c1ccc([Cl,Br,I,F])c([Cl,Br,I,F])c1[Cl,Br,I,F]",
    "c1cc([Cl,Br,I,F])cc([Cl,Br,I,F])c1[Cl,Br,I,F]",
    "[CR1]1[CR1][CR1][CR1][CR1][CR1][CR1]1",
    "[CR1]1[CR1][CR1]cc[CR1][CR1]1",
    "[CR2]1[CR2][CR2][CR2][CR2][CR2][CR2][CR2]1",
    "[CR2]1[CR2][CR2]cc[CR2][CR2][CR2]1",
    "[CH2R2]1N[CH2R2][CH2R2][CH2R2][CH2R2][CH2R2]1",
    "[CH2R2]1N[CH2R2][CH2R2][CH2R2][CH2R2][CH2R2][CH2R2]1",
    "C#C",
    "[OR2,NR2]@[CR2]@[CR2]@[OR2,NR2]@[CR2]@[CR2]@[OR2,NR2]",
    "[$([N+R]),$([n+R]),$([N+]=C)][O-]",
    "[C,c]=N[OH]",
    "[C,c]=NOC=O",
    "[C,c](=O)[CX4,CR0X3,O][C,c](=O)",
    "c1ccc2c(c1)ccc(=O)o2",
    "[O+,o+,S+,s+]",
    "N=C=O",
    "[NX3,NX4][F,Cl,Br,I]",
    "c1ccccc1OC(=O)[#6]",
    "[CR0]=[CR0][CR0]=[CR0]",
    "[C+,c+,C-,c-]",
    "N=[N+]=[N-]",
    "C12C(NC(N1)=O)CSC2",
    "c1c([OH])c([OH,NH2,NH])ccc1",
    "P",
    "[N,O,S]C#N",
    "C=C=O",
    "[Si][F,Cl,Br,I]",
    "[SX2]O",
    "[SiR0,CR0](c1ccccc1)(c2ccccc2)(c3ccccc3)",
    "O1CCCCC1OC2CCC3CCCCC3C2",
    "N=[CR0][N,n,O,S]",
    "[cR2]1[cR2][cR2]([Nv3X3,Nv4X4])[cR2][cR2][cR2]1[cR2]2[cR2][cR2][cR2]([Nv3X3,Nv4X4])[cR2][cR2]2",
    "C=[C!r]C#N",
    "[cR2]1[cR2]c([N+0X3R0,nX3R0])c([N+0X3R0,nX3R0])[cR2][cR2]1",
    "[cR2]1[cR2]c([N+0X3R0,nX3R0])[cR2]c([N+0X3R0,nX3R0])[cR2]1",
    "[cR2]1[cR2]c([N+0X3R0,nX3R0])[cR2][cR2]c1([N+0X3R0,nX3R0])",
    "[OH]c1ccc([OH,NH2,NH])cc1",
    "c1ccccc1OC(=O)O",
    "[SX2H0][N]",
    "c12ccccc1(SC(S)=N2)",
    "c12ccccc1(SC(=S)N2)",
    "c1nnnn1C=O",
    "s1c(S)nnc1NC=O",
    "S1C=CSC1=S",
    "C(=O)Onnn",
    "OS(=O)(=O)C(F)(F)F",
    "N#CC[OH]",
    "N#CC(=O)",
    "S(=O)(=O)C#N",
    "N[CH2]C#N",
    "C1(=O)NCC1",
    "S(=O)(=O)[O-,OH]",
    "NC[F,Cl,Br,I]",
    "C=[C!r]O",
    "[NX2+0]=[O+0]",
    "[OR0,NR0][OR0,NR0]",
    "C(=O)O[C,H1].C(=O)O[C,H1].C(=O)O[C,H1]",
    "[CX2R0][NX3R0]",
    "c1ccccc1[C;!R]=[C;!R]c2ccccc2",
    "[NX3R0,NX4R0,OR0,SX2R0][CX4][NX3R0,NX4R0,OR0,SX2R0]",
    "[s,S,c,C,n,N,o,O]~[n+,N+](~[s,S,c,C,n,N,o,O])(~[s,S,c,C,n,N,o,O])~[s,S,c,C,n,N,o,O]",
    "[s,S,c,C,n,N,o,O]~[nX3+,NX3+](~[s,S,c,C,n,N])~[s,S,c,C,n,N]",
    "[*]=[N+]=[*]",
    "[SX3](=O)[O-,OH]",
    "N#N",
    "F.F.F.F",
    "[R0;D2][R0;D2][R0;D2][R0;D2]",
    "[cR,CR]~C(=O)NC(=O)~[cR,CR]",
    "C=!@CC=[O,S]",
    "[#6,#8,#16][C,c](=O)O[C,c]",
    "c[C;R0](=[O,S])[C,c]",
    "c[SX2][C;!R]",
    "C=C=C",
    "c1nc([F,Cl,Br,I,S])ncc1",
    "c1ncnc([F,Cl,Br,I,S])c1",
    "c1nc(c2c(n1)nc(n2)[F,Cl,Br,I])",
    "[C,c]S(=O)(=O)c1ccc(cc1)F",
    "[15N]",
    "[13C]",
    "[18O]",
    "[34S]",
]

StructuralAlerts = []
for smarts in StructuralAlertSmarts:
    StructuralAlerts.append(Chem.MolFromSmarts(smarts))


# ADS parameters for the 8 molecular properties: [row][column]
#     rows[8]:     MW, ALOGP, HBA, HBD, PSA, ROTB, AROM, ALERTS
#     columns[7]: A, B, C, D, E, F, DMAX
# ALOGP parameters from Gregory Gerebtzoff (2012, Roche)
pads1 = [
    [
        2.817065973,
        392.5754953,
        290.7489764,
        2.419764353,
        49.22325677,
        65.37051707,
        104.9805561,
    ],
    [
        0.486849448,
        186.2293718,
        2.066177165,
        3.902720615,
        1.027025453,
        0.913012565,
        145.4314800,
    ],
    [
        2.948620388,
        160.4605972,
        3.615294657,
        4.435986202,
        0.290141953,
        1.300669958,
        148.7763046,
    ],
    [
        1.618662227,
        1010.051101,
        0.985094388,
        0.000000001,
        0.713820843,
        0.920922555,
        258.1632616,
    ],
    [
        1.876861559,
        125.2232657,
        62.90773554,
        87.83366614,
        12.01999824,
        28.51324732,
        104.5686167,
    ],
    [
        0.010000000,
        272.4121427,
        2.558379970,
        1.565547684,
        1.271567166,
        2.758063707,
        105.4420403,
    ],
    [
        3.217788970,
        957.7374108,
        2.274627939,
        0.000000001,
        1.317690384,
        0.375760881,
        312.3372610,
    ],
    [
        0.010000000,
        1199.094025,
        -0.09002883,
        0.000000001,
        0.185904477,
        0.875193782,
        417.7253140,
    ],
]
# ALOGP parameters from the original publication
pads2 = [
    [
        2.817065973,
        392.5754953,
        290.7489764,
        2.419764353,
        49.22325677,
        65.37051707,
        104.9805561,
    ],
    [
        3.172690585,
        137.8624751,
        2.534937431,
        4.581497897,
        0.822739154,
        0.576295591,
        131.3186604,
    ],
    [
        2.948620388,
        160.4605972,
        3.615294657,
        4.435986202,
        0.290141953,
        1.300669958,
        148.7763046,
    ],
    [
        1.618662227,
        1010.051101,
        0.985094388,
        0.000000001,
        0.713820843,
        0.920922555,
        258.1632616,
    ],
    [
        1.876861559,
        125.2232657,
        62.90773554,
        87.83366614,
        12.01999824,
        28.51324732,
        104.5686167,
    ],
    [
        0.010000000,
        272.4121427,
        2.558379970,
        1.565547684,
        1.271567166,
        2.758063707,
        105.4420403,
    ],
    [
        3.217788970,
        957.7374108,
        2.274627939,
        0.000000001,
        1.317690384,
        0.375760881,
        312.3372610,
    ],
    [
        0.010000000,
        1199.094025,
        -0.09002883,
        0.000000001,
        0.185904477,
        0.875193782,
        417.7253140,
    ],
]


def ads(x, a, b, c, d, e, f, dmax):
    return (
        a
        + (
            b
            / (1 + exp(-1 * (x - c + d / 2) / e))
            * (1 - 1 / (1 + exp(-1 * (x - c - d / 2) / f)))
        )
    ) / dmax


def properties(mol):
    """
    Calculates the properties that are required to calculate the QED descriptor.
    """
    matches = []
    if mol is None:
        raise WrongArgument("properties(mol)", "mol argument is 'None'")
    x = [0] * 9
    x[0] = Descriptors.MolWt(mol)  # MW
    x[1] = Descriptors.MolLogP(mol)  # ALOGP
    for hba in Acceptors:  # HBA
        if mol.HasSubstructMatch(hba):
            matches = mol.GetSubstructMatches(hba)
            x[2] += len(matches)
    x[3] = Descriptors.NumHDonors(mol)  # HBD
    x[4] = Descriptors.TPSA(mol)  # PSA
    x[5] = Descriptors.NumRotatableBonds(mol)  # ROTB
    x[6] = Chem.GetSSSR(Chem.DeleteSubstructs(deepcopy(mol), AliphaticRings))  # AROM
    for alert in StructuralAlerts:  # ALERTS
        if mol.HasSubstructMatch(alert):
            x[7] += 1
    ro5_failed = 0
    if x[3] > 5:
        ro5_failed += 1  # HBD
    if x[2] > 10:
        ro5_failed += 1  # HBA
    if x[0] >= 500:
        ro5_failed += 1
    if x[1] > 5:
        ro5_failed += 1
    x[8] = ro5_failed
    return x


def qed(w, p, gerebtzoff):
    d = [0.00] * 8
    if gerebtzoff:
        for i in range(0, 8):
            d[i] = ads(
                p[i],
                pads1[i][0],
                pads1[i][1],
                pads1[i][2],
                pads1[i][3],
                pads1[i][4],
                pads1[i][5],
                pads1[i][6],
            )
    else:
        for i in range(0, 8):
            d[i] = ads(
                p[i],
                pads2[i][0],
                pads2[i][1],
                pads2[i][2],
                pads2[i][3],
                pads2[i][4],
                pads2[i][5],
                pads2[i][6],
            )
    t = 0.0
    for i in range(0, 8):
        t += w[i] * log(d[i])
    return exp(t / sum(w))


def weights_max(mol, gerebtzoff=True, props=False):
    """
    Calculates the QED descriptor using maximal descriptor weights.
    If props is specified we skip the calculation step and use the props-list of properties.
    """
    if not props:
        props = properties(mol)
    return qed([0.50, 0.25, 0.00, 0.50, 0.00, 0.50, 0.25, 1.00], props, gerebtzoff)


def weights_mean(mol, gerebtzoff=True, props=False):
    """
    Calculates the QED descriptor using average descriptor weights.
    If props is specified we skip the calculation step and use the props-list of properties.
    """
    if not props:
        props = properties(mol)
    return qed([0.66, 0.46, 0.05, 0.61, 0.06, 0.65, 0.48, 0.95], props, gerebtzoff)


def weights_none(mol, gerebtzoff=True, props=False):
    """
    Calculates the QED descriptor using unit weights.
    If props is specified we skip the calculation step and use the props-list of properties.
    """
    if not props:
        props = properties(mol)
    return qed([1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00], props, gerebtzoff)


def default(mol, gerebtzoff=True):
    """
    Calculates the QED descriptor using average descriptor weights and Gregory Gerebtzoff parameters.
    """
    return weights_mean(mol, gerebtzoff)


if __name__ == "__main__":
    parser = argparse.ArgumentParser()
    parser.add_argument(
        "-i", "--input", required=True, help="path to the input file name"
    )
    parser.add_argument(
        "-m",
        "--method",
        dest="method",
        choices=["max", "mean", "unweighted"],
        default="mean",
        help="Specify the method you want to use.",
    )
    parser.add_argument(
        "--iformat", help="Input format. It must be supported by openbabel."
    )
    parser.add_argument(
        "-o",
        "--outfile",
        type=argparse.FileType("w+"),
        default=sys.stdout,
        help="path to the result file, default it sdtout",
    )
    parser.add_argument(
        "--header",
        dest="header",
        action="store_true",
        default=False,
        help="Write header line.",
    )

    args = parser.parse_args()

    # Elucidate filetype and open supplier
    ifile = os.path.abspath(args.input)
    if not os.path.isfile(ifile):
        print("Error: ", ifile, " is not a file or cannot be found.")
        sys.exit(1)
    if not os.path.exists(ifile):
        print("Error: ", ifile, " does not exist or cannot be found.")
        sys.exit(1)
    if not os.access(ifile, os.R_OK):
        print("Error: ", ifile, " is not readable.")
        sys.exit(1)

    if not args.iformat:
        # try to guess the filetype
        filetype = check_filetype(ifile)
    else:
        filetype = args.iformat  # sdf or smi

    """
        We want to store the original SMILES in the output. So in case of a SMILES file iterate over the file and convert each line separate.
    """
    if filetype == "sdf":
        supplier = Chem.SDMolSupplier(ifile)
        # Process file
        if args.header:
            args.outfile.write(
                "MW\tALOGP\tHBA\tHBD\tPSA\tROTB\tAROM\tALERTS\tLRo5\tQED\tNAME\n"
            )
        count = 0
        for mol in supplier:
            count += 1
            if mol is None:
                print(
                    "Warning: skipping molecule ", count, " and continuing with next."
                )
                continue
            props = properties(mol)

            if args.method == "max":
                calc_qed = weights_max(mol, True, props)
            elif args.method == "unweighted":
                calc_qed = weights_none(mol, True, props)
            else:
                calc_qed = weights_mean(mol, True, props)

            args.outfile.write(
                "%.2f\t%.3f\t%d\t%d\t%.2f\t%d\t%d\t%d\t%s\t%.3f\t%-s\n"
                % (
                    props[0],
                    props[1],
                    props[2],
                    props[3],
                    props[4],
                    props[5],
                    props[6],
                    props[7],
                    props[8],
                    calc_qed,
                    mol.GetProp("_Name"),
                )
            )
    elif filetype == "smi":
        supplier = Chem.SmilesMolSupplier(ifile, " \t", 0, 1, False, True)

        # Process file
        if args.header:
            args.outfile.write(
                "MW\tALOGP\tHBA\tHBD\tPSA\tROTB\tAROM\tALERTS\tLRo5\tQED\tNAME\tSMILES\n"
            )
        count = 0
        for line in open(ifile):
            tokens = line.strip().split("\t")
            if len(tokens) > 1:
                smiles, title = tokens
            else:
                smiles = tokens[0]
                title = ""
            mol = Chem.MolFromSmiles(smiles)
            count += 1
            if mol is None:
                print(
                    "Warning: skipping molecule ", count, " and continuing with next."
                )
                continue
            props = properties(mol)

            if args.method == "max":
                calc_qed = weights_max(mol, True, props)
            elif args.method == "unweighted":
                calc_qed = weights_none(mol, True, props)
            else:
                calc_qed = weights_mean(mol, True, props)

            args.outfile.write(
                "%.2f\t%.3f\t%d\t%d\t%.2f\t%d\t%d\t%d\t%s\t%.3f\t%-s\t%s\n"
                % (
                    props[0],
                    props[1],
                    props[2],
                    props[3],
                    props[4],
                    props[5],
                    props[6],
                    props[7],
                    props[8],
                    calc_qed,
                    title,
                    smiles,
                )
            )
    else:
        sys.exit("Error: unknown file-type: %s" % filetype)