B 3Rca@sdZddlmZddlmZmZddlmZddlmZddl m Z ddl m Z ddl mZmZdd lmZdd lZdd lZdd lZdd lmZdd lZdd lZdd lZdd lmZed ZeddddgZedZddZGdddZ GdddZ!e"Z#ej$e#_%de#_&ej'j(e#_)de#_*e#fe!dddZ+d+ddZ,d,e!eed d fdd d!Z-d"d#Z.e!d$d%d&Z/d'd(Z0d)d*Z1d S)-u Free Wilson Analysis A Mathematical Contribution to Structurre-Activity Studies Spencer M, Free, Jr and James W Wilson Journal of Medicinal Chemistry Vol 7, Number 4, July 6, 1964 Basic usage: get a scaffold (or scaffolds) some compounds and their scores, then run freewilson: Read some example data: ``` >>> import os, sys >>> DATA_PATH = "data" >>> smilesfile = os.path.join(DATA_PATH, "CHEMBL2321810.smi") >>> scaffoldfile = os.path.join(DATA_PATH, "CHEMBL2321810_scaffold.mol") >>> csvfile = os.path.join(DATA_PATH, "CHEMBL2321810_act.csv") >>> mols = [] >>> for line in open(smilesfile): ... smiles, name = line.strip().split() ... m = Chem.MolFromSmiles(smiles) ... m.SetProp("_Name", name) ... mols.append(m) >>> scaffold = Chem.MolFromMolBlock(open(scaffoldfile).read()) >>> data = {k:float(v) for k,v in list(csv.reader(open(csvfile)))[1:]} >>> scores = [data[m.GetProp("_Name")] for m in mols] ``` And do the decomposition: ``` >>> from freewilson import FWDecompose, FWBuild, predictions_to_csv >>> decomp = FWDecompose(scaffold, mols, scores) ``` Scores need to be in an appropriate form for regression analysis, i.e. pIC50s as opposed to IC50s. Enumerations are a lot faster if you know a prediction cutoff or molecular weight cutoff. Scaffolds can be any scaffold or smarts pattern or list of either. The system automatically enumerates matching cores in the analysis. To see if the decomposition is useful, check the r squared value ``` >>> print(f"Training R^2 is {decomp.r2:0.2f}") Training R^2 is 0.81 ``` Finally you can build the decomposition into new molecules: ``` >>> for pred in FWBuild(decomp): # doctest: +SKIP ... print(pred.smiles, pred.prediction) ``` Now this builds both well and poorly predicted molecules. To prevent this you can use the following filters while building: 1. pred_filter: given a prediction, return True to keep the molecule 2. hvy_filter: given a heavy atom count, return True to keep the molecule 3. mw_filter: given a molecular weight, return True to keep the molecule 4. mol_filter: given a sanitized molecule, return True to keep the molecule Here are some examples (see using Molecular Filters below) ``` >>> preds = FWBuild(decomp, ... pred_filter=lambda x: x > 8, ... mw_filter=lambda mw: 100>> predictions_to_csv(sys.stdout, decomp, preds) ``` More Info --------- You can also get some more information by setting logging to INFO ``` >>> import logging >>> logging.getLogger().setLevel(logging.INFO) >>> preds = list(FWBuild(decomp, pred_filter=lambda x: x > 8)) ``` You'll see something like this: ``` INFO:root:Enumerating rgroups with no broken cycles... INFO:root: Core 1 INFO:root: R1 73 INFO:root: R10 2 INFO:root: R3 445 100%|███████████████████████████████████████████████████████| 64970/64970 [00:05<00:00, 11247.38it/s] INFO:root:Wrote 3 results out of 64970 In Training set: 628 Bad MW: 0 Bad Pred: 64339 Bad Filters: 0 Bad smi: 0 min mw: 380.429 max mw: 772.4030000000001 min pred: 2.8927324757327666 max pred: 8.148182660251193 ``` Using FMCS to find a scaffold ----------------------------- If you don't have a good feel for the dataset, you can generate the scaffold by using the rdFMCS package. The following tries to find a decent MCS that covers 80% of the input structures ``` >>> from rdkit.Chem import rdFMCS >>> mcs = rdFMCS.FindMCS(mols[0:8], threshold=0.8, atomCompare=rdFMCS.AtomCompare.CompareAny, ... completeRingsOnly=True) # doctest: +SKIP >>> decomp = FWDecompose(mcs.queryMol, mols, scores) # doctest: +SKIP ``` )rdRGroupDecomposition)molzip Descriptors)rdBase)Chem)Ridge)r2_score) defaultdict namedtuple)tqdmN)List) GeneratorZ freewilsonFreeWilsonPrediction predictionsmilesrgroupsz \*:([0-9]+)cCst}g}x@|dD]2}|ddkr>||kr4qn ||||qWd|}ttj|dd}t |}|d}i}t t }x6|D].} | |  } | r| || <|| d7<qWt|d} g} xz|ddD]j} xd| D]\} | |  } | r|| d}|| <|dkr| | | | || | f| d7} qWqWx| D]~\}} | |}t|}t|dkrR|d|}td}||} | | }|| |||tjjqRWt|S) z~Fix a rgroup smiles for molzip, note that the core MUST come first in the smiles string, ala core.rgroup1.rgroup2 ... .*F)sanitizerN)setsplitcountaddappendjoinrZRWMol MolFromSmilesZ GetMolFragsr intZGetAtomWithIdxZ GetAtomMapNummaxZ SetAtomMapNumlistZGetBondslenZ GetOtherAtomZAtomZAddAtomZAddBondZGetIdxZBondTypeZSINGLEr)rZdupesZslsmZfragsZcoreZatommapsZcountsidxZatommapZ next_atommapZ add_atommapZfragmentrZatomidxatomZbondsZoatomZxatomr&,share/RDKit/Contrib/FreeWilson/freewilson.py molzip_smisR             r(c@s*eZdZdZd ddZddZddZdS) RGroupaFreeWilson RGroup smiles - isomeric smiles of the rgroup rgroup - Rgroup name ['Core', 'R1', 'R2'...] count - number of molecules with the rgroup coefficient - ridge regression coefficient idx - one-hot encoding for the rgroup Nc Cs||_||_||_||_||_tddtt|D|_ t |}|rdt ||_t ||_nt }t j |dd}|rd|_d|_xB|D]6}|} |j|| 7_| dkr|jd7_qWt ||_n d|_d|_dS)NcSsg|] }t|qSr&)r).0xr&r&r' sz#RGroup.__init__..F)rgrr)rrgroupr coefficientr$tuplesorteddummypatfindalldummiesrrrMolWtmwHeavyAtomCounthvyctZGetPeriodicTableZGetAtomsZ GetAtomicNumZGetAtomicWeight) selfrr-rr.r$r#tabler%Z atomicnumr&r&r'__init__s.  zRGroup.__init__c CsBdt|jdt|jdt|jdt|jdt|jd S)NzRGroup(smiles=z , rgroup=z, count=z, coefficient=z, idx=))reprrr-rr.r$)r8r&r&r'__str__szRGroup.__str__cCs|S)N)r=)r8r&r&r'__repr__szRGroup.__repr__)N)__name__ __module__ __qualname____doc__r:r=r>r&r&r&r'r)s r)c@seZdZdZddZdS)FreeWilsonDecompositionaFreeWilson decomposition rgroups - dictionary of rgroup to list of RGroups i.e. {'Core': [RGroup(...), ...] 'R1': [ RGroup(...), RGroup(...)], } rgroup_to_descriptor_idx - one hot encoding of smiles to descriptor index fitter - scikit learn compatible fitter used for regression N - number of rgroups r2 - regression r squared descriptors - set of the descriptors for molecules in the training set used to not enumerate existing molecules row_decomposition - original rgroup decomposition (With row key 'molecule' is an rdkit molecule) c CsL||_||_||_t||_||_tdd|D|_||_||_ ||_ dS)NcSsg|] }t|qSr&)r/)r*dr&r&r'r,sz4FreeWilsonDecomposition.__init__..) rrgroup_to_descriptor_idxfitterr!Nr2r descriptorsrow_decomposition num_trainingnum_reconstructed) r8rrErFrHrIrJrKrLr&r&r'r:s z FreeWilsonDecomposition.__init__N)r?r@rArBr:r&r&r&r'rCs rCF)returnc Csg}g}i}tt}t|t|kr@tdt|dt|tdt|dt||}g} g} xNtt t ||D]8\} \} } | | dkr~| | | | | | q~W| tdt|dt||stdd S|jd d }td t}tt}d}xtt || D]\}\}}g}x^|D]R\}}| |||d 7<||krFt|||<|| t||ddqFW||d<d|}yt}t|} |d 7}Wn"tdt| ||YnXq,Wtdt|td|d||dkr6td||xbt | |D]T\} }| |d<dgt|}| |x(|D]}||krtd |||<qtWqBWt|t|kstdt|dt|dtdt}|||| |}t!||}td|td|j"dxL|D]@}x8|D]0}||j#|_$|j%||j#|_&||j#|_'q2Wq(Wt(||||||||S)ad Perform a free wilson analysis : param scaffolds : scaffold or list of scaffolds to use for the rgroup decomposition : param mols : molecules to decompose : param scores : list of floating point numbers for the regression ( you may need convert these to their logs in some cases) : param decomp_params : RgroupDecompositionParams default [ default_decomp_params = rdkit.Chem.rdRGroupDecomposition.RGroupDecompositionParameters() default_decomp_params.matchingStrategy = rgd.GA default_decomp_params.onlyMatchAtRGroups = False ] If you only want to decompose on specific group locations set onlyMatchAtRGroups to True >>> from rdkit import Chem >>> from freewilson import FWBuild, FWDecompose >>> from rdkit.Chem import Descriptors >>> scaffold = Chem.MolFromSmiles("c1cccnc1") >>> mols = [Chem.MolFromSmiles("c1cccnc1"+"C"*(i+1)) for i in range(100)] >>> scores = [Descriptors.MolLogP(m) for m in mols] >>> fw = FWDecompose(scaffold, mols, scores) >>> for pred in FWBuild(fw): ... print(pred) For an easy way to report predictions see >>> from freewilson import FWBuild, predictions_to_csv >>> import sys >>> predictions_to_csv(sys.stdout, FWBuild(fw)) See FWBuild docs to see how to filter predictions, molecular weight or molecular properties. z>The number of molecules must match the number of scores #mols z #scores z Decomposing z molecules...rzMatched z out of z(No scaffolds matched the input moleculesNT)ZasSmileszGet unique rgroups...rZ original_idxrzfailed:zDescriptor size zReconstructed z6Could only reconstruct %s out of %s training moleculesZmoleculezNumber of descriptors(z)) doesn't match number of matcved scores(r;zRidge Regressing...zR2 z Intercept = z.2f))r r r! ValueErrorloggerinforgdZRGroupDecomposition enumerater zipZAddrZProcesserrorZGetRGroupsAsRowsr BlockLogsritemsr)rr(printr MolToSmilesloggingZwarningvaluesAssertionErrorrZfitpredictrZ intercept_rrZcoef_r.r$rC)Z scaffoldsmolsscoresZ decomp_paramsrIZmatched_scoresZ rgroup_idxrZ decomposerZmatchedZmatched_indicesimolZscore decompositionblockerZ rgroup_countsrLZnum_molsrowr$Z row_smilesr-rZ reconstructedZ descriptorZlmZpredsrHZ sidechainsr&r&r' FWDecompose-s#            "   0      rdc%(cs|j}|j}d}x|D]} |t| 9}qWd} d} d} d} d}d}d}d}d}d}d}d}d}d}|dpld}xtttj||dD]\}}|r||dkrtjd| d |d | d | d | d |d|d|d|d|d|t j dd}d}dg|}x,|D]$}d||j <||j 7}||j 7}q Wt||jkrN| d7} qt||}t||}|r|||s|| d7} qt||}t||}|r||s|d7}q||gd}t||}t||}|r||s| d7} q|d7}tdd|D} dtdd|D}!y t|!}"Wn|d7}wYnXd}#|r\||"s\|d7}qt|"}$t||$|V| d7} qWtd| d |d| d | d | d |d|d|d|d|d | d |d|d|d|d|d|d|d|&dS)Nrrig _g _B )ZtotalzWrote z results out of z In Training Setz Bad MW: z Bad Pred: z Bad Filters: z Bad smi: z min mw: z max mw: z min pred: z max pred: )filecSsg|] }|jqSr&)r)r*gr&r&r'r,sz_enumerate..rcSsg|] }|jqSr&)r)r*rgr&r&r'r,sFz In Training set: z Bad HVY: z min hvy: z max hvy: )rGrFr!r rR itertoolsproductrYdebugsysstderrr$r5r7r/rIminrr\rrr(rrXrrP)%rfw mw_filter hvy_filter pred_filter mol_filterrGrFZ num_productsrZwroteZin_training_setZ rejected_predZ rejected_mwZ rejected_hvyZ good_predZrejected_filtersZ rejected_badZmin_mwZmax_mwZmin_hvyZmax_hvyZmax_predZmin_predZdeltar_groupsr5hvyrIrgpredrZsmir`ZrejectedZout_smir&r&r' _enumerates  $R             rw)rnrMc cs>t}t}tt}tt}x|jD]r\} } | dkrD| || <q*|| } xN| D]F} | | qRt| j dkr| | j || j | qR| | qRWq*Wt dx.|D]"\} }t d| dt|qWddt |D}x$t||||||dD] }|VqWt}x4|jD]*} | dd kr$| t| dd q$W|rbt d x|D]}|t|}d|jdi}|||d d dd|D<x"|D]}d |} || || <qWx0|D]$\} }t d| dt|qWx6tddt |D|||||dD] }|Vq$WqhWd S)aEnumerate the freewilson decomposition and return their predictions :param fw: FreeWilsonDecomposition generated from FWDecompose :param pred_filter: return True if the prediction is in a desireable range e.g. lambda pic50: pic50>=8 :param mw_filter: return True if the enumerated molecular weight is in a desireable rrange e.g. lambda mw: 150 < mw < 550 :param hvy_filter: return True if the enumerated heavy couont is in a desireable rrange e.g. lambda hvy: 10 < hvy < 50 :param mol_filter: return True if the molecule is ok to be enumerated e.g. lambda mol: -3 < Descriptors.MolLogp(mol) < 5 Corerz,Enumerating rgroups with no broken cycles... cSsg|] \}}|qSr&r&)r*keyr-r&r&r'r,szFWBuild..)rqrorprrrRNz)Enumerating rgroups with broken cycles...zR%srcSsg|] }t|qSr&)str)r*r+r&r&r'r,5scSsg|] \}}|qSr&r&)r*rzr-r&r&r'r,=s)rrUrr r rrVrr!r3rrYrPr0rwrr)rnrqrorprrrbZcyclesZrgroups_no_cyclesZ rgroup_cyclesrzr-Z no_cyclesrgkvrZresindicesZrgroup_indicesZmissingr#r&r&r'FWBuildsb              rcCs(|ddkr dt|ddfS|dfS)z-Sort groups like R1 R2 R10 not R1 R10 R2 rr{rN)r)rsr&r&r' _rgroup_sortEs r)rac Csd}x|D]}|st}x|jD]}||q Wt|td}i}xt|D]\}} ||| <qLWt|}ddgddt|D} | | dgt |} x|jD]} | j | || j <qW|j t |jg| } | | q W| S)zOutput predictions in csv format to the output stream :param outstream: output stream to write results :param decomposition: freewillson decomposition :param predictions: list of Predictions to output N)rzrrcSsg|]}|dqS)Z_smilesr&)r*rgr&r&r'r,^sz&predictions_to_csv..)rrrr0rrRcsvwriterr Zwriterowr!rr-r<r) Z outstreamraZ predictionsrrvrr-lookupr_rheaderr"rcr&r&r'predictions_to_csvKs&       rcs"ddlmddlmtdddgks.ttdddgksDttdddgksZttd ddgksptd }fd d td D}fdd |D}t |||}xnt |ddddddfdddD]D}t }x|j D]}| |jqWt|}t|dddgkstqWdS)Nr)r)rzC[*:1]N.[H][*:2]12zC[*:1]N.[HH][*:2]zC[*:1]N.[2H][*:2]zC[*:1]N.[CH2][*:2]z[*:2]c1cccnc1[*:1]cs.g|]&}d|ddd|dqS)rGrZc1cccnc1C)r)r*r_)rr&r'r,rsz#test_freewilson..recsg|]}|qSr&)MolLogP)r*r#)rr&r'r,sscSsd|kodkSS)Nr&)r+r&r&r'wz!test_freewilson..cSsd|kodkSS)Ndir&)r5r&r&r'rxrcSsd|kodkSS)Nre2r&)rur&r&r'ryrcsd|kodkSS)Nrr)r)r#)rr&r'rzr)rqrorprrrxZR1ZR2)rdkitr rdkit.Chemrr1r2r[rrangerdrrrrr-r0r )Zscaffoldr]r^rnrvrZ sidechainr&)rrr'test_freewilsonhs(     rcCsd}t|dddd}t|}|jt|ks2t|jt|ksFt|j dksTtd}t|dddd}|jdksvt|jd kst|j dkstd }t|dddd}|jd kst|jd kst|j t kstdS) Nz[*:2]c1ccccc1[*:1]rxrg?)rr-rr.)rrz[*:2]cccccc[*:1]gMbX9R@ZNoper) r)rrr5rr4r[r7r6r3r/)rrr`r&r&r' test_rgroupss  r)NNNN)NNNN)2rBrrrQrrrrrZsklearn.linear_modelrZsklearn.metricsr collectionsr r r rhZmathrktypingr rYrrer Z getLoggerrOrcompiler1r(r)rCZRGroupDecompositionParametersZdefault_decomp_paramsZGAZmatchingStrategyZonlyMatchAtRGroupsZ RGroupScoreZFingerprintVarianceZ scoreMethodZremoveHydrogensPostMatchrdrwrrrrrr&r&r&r'sJ          2- { MM