MolFilterGAN/ScoreFunc.py at main · xiaohongniua/MolFilterGAN · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
# -*- coding: utf-8 -*-
"""
Created on Wed Nov  6 21:01:43 2019

@author: xhliu
"""
try:
    import rdkit
    from rdkit import Chem
    from rdkit import rdBase
    from rdkit.Chem.rdchem import HybridizationType
    from rdkit.Chem import rdMolDescriptors as ds
    import numpy as np
    from rdkit import RDConfig
    from rdkit.Chem import ChemicalFeatures,Lipinski
    from rdkit.Chem.rdchem import BondType as BT
    import sascorer
    from rdkit.Chem import ChemicalFeatures,Descriptors,rdMolDescriptors,QED
    rdBase.DisableLog('rdApp.error')
except ImportError:
    rdkit = None


def rule5(infile,outfile):
    a=open(infile)
    b=open(outfile,'w')
    for i in a:
        smiles=i.rstrip('\n')
        mol=Chem.MolFromSmiles(smiles)
        if mol:
            logp=Descriptors.MolLogP(mol)
            mw=Descriptors.MolWt(mol)
            hbd=rdMolDescriptors.CalcNumLipinskiHBD(mol)
            hda=rdMolDescriptors.CalcNumLipinskiHBA(mol)
            rb=rdMolDescriptors.CalcNumRotatableBonds(mol)
            b.write(smiles+','+str(logp)+','+str(mw)+','+str(hbd)+','+str(hda)+','+str(rb)+'\n')
    a.close()
    b.close()


def ca_sa(insmi):
    mol = Chem.MolFromSmiles(insmi)
    s = sascorer.calculateScore(mol)
    return s


def ca_qed(insmi):
    mol = Chem.MolFromSmiles(insmi)
    q = QED.qed(mol)
    return q


def CalculateZagreb1(mol):
    """
    #################################################################
    Calculation of Zagreb index with order 1 in a molecule

    ---->ZM1

    Usage:

        result=CalculateZagreb1(mol)

        Input: mol is a molecule object

        Output: result is a numeric value
    #################################################################
    """

    deltas = [x.GetDegree() for x in mol.GetAtoms()]
    return sum(np.array(deltas) ** 2)


def CalculateQuadratic(mol):
    """
    #################################################################
    Calculation of Quadratic index in a molecule

    ---->Qindex

    Usage:

        result=CalculateQuadratic(mol)

        Input: mol is a molecule object

        Output: result is a numeric value
    #################################################################
    """
    deltas = [x.GetDegree() for x in mol.GetAtoms()]
    M = sum(np.array(deltas) ** 2)
    N = mol.GetNumAtoms()
    return 3 - 2 * N + M / 2.0


def ca_MCE18(in_smi):
    mce_score = 0
    mol = Chem.MolFromSmiles(in_smi)
    AR = ds.CalcNumAromaticRings(mol)
    # print('AR:%s'%AR)
    if AR > 0:
        mce_score += 1
    NAR = ds.CalcNumAliphaticRings(mol)
    # print('NAR:%s'%NAR)
    if NAR > 0:
        mce_score += 1
    chiral = Chem.FindMolChiralCenters(mol)
    # print('chiral:%s' % chiral)
    if len(chiral) > 0:
        mce_score += 1
    spiro = ds.CalcNumSpiroAtoms(mol)
    # print('spiro:%s' % spiro)
    if spiro > 0:
        mce_score += 1

    sp3 = ds.CalcFractionCSP3(mol)
    # print('csp3:%s' % sp3)
    atom_num = 0
    C_num = 0
    csp3_num = 0
    csp3_cyc_num = 0
    csp3_acyc_num = 0
    for atom in mol.GetAtoms():
        atom_num += 1
        if atom.GetSymbol() == 'C':
            C_num += 1
            if atom.GetHybridization() == Chem.rdchem.HybridizationType.SP3:
                csp3_num += 1
                if atom.IsInRing():
                    csp3_cyc_num += 1
                else:
                    csp3_acyc_num += 1
            else:
                pass
        else:
            pass
    # csp3_C_num = csp3_cyc_num + csp3_acyc_num
    # cyc = csp3_cyc_num/csp3_num
    # acyc = csp3_acyc_num/csp3_num
    cyc = csp3_cyc_num / C_num
    acyc = csp3_acyc_num/C_num
    NCSPTR = (sp3 + cyc - acyc)/(1+sp3)
    mce_score += NCSPTR
    # print(mce_score)
    normed_q = CalculateQuadratic(mol)
    mce_score = mce_score * normed_q
    return mce_score


def ca_fsp3(insmi):
    mol = Chem.MolFromSmiles(insmi)
    p3 = Lipinski.FractionCSP3(mol)
    return p3