import logging import csv import re import platform import io from subprocess import Popen, PIPE #from openbabel import openbabel import pdb from rdkit.Chem import rdchem if platform.system() == 'Windows': CXCALC_BIN = 'C:\\Users\\vuu10\\AppData\\Local\\Programs\\ChemAxon\\MarvinSuite\\bin\\cxcalc.exe' #CXCALC_BIN = 'C:\\Program Files (x86)\\ChemAxon\\MarvinBeans\\bin\\cxcalc.bat' use_shell_for_echo = True else: CXCALC_BIN = 'cxcalc' use_shell_for_echo = False MID_PH = 7.0 N_PKAS = 20 class ChemAxonError(Exception): pass def RunCxcalc(molstring, args): # pdb.set_trace() # with open(platform.DEV_NULL, 'w') as dev_null: try: logging.debug("INPUT: echo %s | %s" % (molstring, ' '.join([CXCALC_BIN] + args))) p1 = Popen(["echo", molstring], stdout=PIPE, shell=use_shell_for_echo) # p2 = Popen([CXCALC_BIN] + args, stdin=p1.stdout, # executable=CXCALC_BIN, stdout=PIPE, stderr=dev_null, shell=False) p2 = Popen([CXCALC_BIN] + args, stdin=p1.stdout, executable=CXCALC_BIN, stdout=PIPE, shell=False) # p.wait() # os.remove(temp_fname) res = p2.communicate()[0] if p2.returncode != 0: raise ChemAxonError(str(args)) logging.debug("OUTPUT: %s" % res) res = res.decode('utf-8') return res except OSError: raise Exception( "Marvin (by ChemAxon) must be installed to calculate pKa data.") def ParsePkaOutput(s, n_acidic, n_basic): """ Returns: A dictionary that maps the atom index to a list of pKas that are assigned to that atom. """ # s = s.decode('utf-8') atom2pKa = {} pkaline = s.split('\n')[1] splitline = pkaline.split('\t') splitline.pop(0) if n_acidic + n_basic > 0: if len(splitline) != (n_acidic + n_basic + 2): raise ChemAxonError('ChemAxon failed to find any pKas') pKa_list = [] acid_or_base_list = [] for i in range(n_acidic + n_basic): x = splitline.pop(0) if x == '': continue pKa_list.append(float(x)) if i < n_acidic: acid_or_base_list.append('acid') else: acid_or_base_list.append('base') atom_list = splitline.pop(0) if atom_list: # a comma separated list of the deprotonated atoms atom_numbers = [int(y)-1 for y in atom_list.split(',')] for i, j in enumerate(atom_numbers): atom2pKa.setdefault(j, []) atom2pKa[j].append((pKa_list[i], acid_or_base_list[i])) smiles_list = splitline return atom2pKa, smiles_list def GetDissociationConstants_val(molstring, n_acidic=N_PKAS, n_basic=N_PKAS, pH=MID_PH): """ Returns: A pair of (pKa list, major pseudoisomer) - the pKa list is of the pKa values in ascending order. - the major pseudoisomer is a SMILES string of the major species at the given pH. """ args = [] if n_acidic + n_basic > 0: args += ['pka', '-a', str(n_acidic), '-b', str(n_basic), 'majorms', '-M', 'true', '--pH', str(pH)] output = RunCxcalc(molstring, args) atom2pKa, smiles_list = ParsePkaOutput(output, n_acidic, n_basic) all_pKas = [] for pKa_list in list(atom2pKa.values()): all_pKas += [pKa for pKa, _ in pKa_list] return sorted(all_pKas), smiles_list def GetDissociationConstants(molstring, n_acidic=N_PKAS, n_basic=N_PKAS, pH=MID_PH): """ Arguments: molstring - a text description of the molecule (SMILES or InChI) n_acidic - the max no. of acidic pKas to calculate n_basic - the max no. of basic pKas to calculate pH - the pH for which the major pseudoisomer is calculated Returns a pair: (all_pKas, major_ms) - all_pKas is a list of floats (pKa values) - major_ms is a SMILES string of the major pseudoisomer at pH_mid """ all_pKas, smiles_list = GetDissociationConstants_val(molstring, n_acidic, n_basic, pH) major_ms = smiles_list[0] return all_pKas, major_ms def GetFormulaAndCharge(molstring): """ Arguments: molstring - a text description of the molecule (SMILES or InChI) Returns: chemical formula of the molecule """ args = ['formula', 'formalcharge'] output = RunCxcalc(molstring, args) # the output is a tab separated table whose columns are: # id, Formula, Formal charge f = io.StringIO(output) tsv_output = csv.reader(f, delimiter='\t') headers = next(tsv_output) if headers != ['id', 'Formula', 'Formal charge']: raise ChemAxonError( 'cannot get the formula and charge for: ' + molstring) _, formula, formal_charge = next(tsv_output) try: formal_charge = int(formal_charge) except ValueError: formal_charge = 0 return formula, formal_charge def GetAtomBagAndCharge(molstring): formula, formal_charge = GetFormulaAndCharge(molstring) periodic_table = rdchem.GetPeriodicTable() atom_bag = {} for mol_formula_times in formula.split('.'): for times, mol_formula in re.findall('^(\d+)?(\w+)', mol_formula_times): if not times: times = 1 else: times = int(times) for atom, count in re.findall("([A-Z][a-z]*)([0-9]*)", mol_formula): if count == '': count = 1 else: count = int(count) atom_bag[atom] = atom_bag.get(atom, 0) + count * times n_protons = sum([c * periodic_table.GetAtomicNumber(str(elem)) for (elem, c) in atom_bag.items()]) atom_bag['e-'] = n_protons - formal_charge return atom_bag, formal_charge if __name__ == "__main__": logging.getLogger().setLevel(logging.WARNING) from molecule import Molecule compound_list = [ ('D-Erythrulose', 'InChI=1S/C4H8O4/c5-1-3(7)4(8)2-6/h3,5-7H,1-2H2/t3-/m1/s1')] for name, inchi in compound_list: print("Formula: %s\nCharge: %d" % GetFormulaAndCharge(inchi)) diss_table, major_ms = GetDissociationConstants(inchi) m = Molecule.FromSmiles(major_ms) print("Name: %s\nInChI: %s\npKas: %s" % (name, m.ToInChI(), str(diss_table)))