Ignota DSG Notes

# Ignota DSG Team Notes - Day 1 Day 2 notes: https://hackmd.io/0qNBkU0HRXKf3Jrgsl7HiQ ## Outline (thanks Pareena for making the outline!) 1) Introductions – come up with team name for breaking the ice 2) Admin a) Agreeing on a way to meet every day b) Talking about Git, make sure everyone knows how to use it and encourage to push and pull often c) Talk about writing up, making sure everyone knows how to use overleaf d) Maybe make draft slides simultaneously while writing up. e) A group review at the end of the day/beginning of the next day about report and progress 3) Group discussion about challenge a) what people think would work well with the problem b) and what they are keen to work on 4) Agreeing on a broad approach a) Split main problem into smaller chunks b) Talking about the common code we need (way to load and cross fold the data, evaluation of models) 5) Split team into small groeups to work on different things ## Organisational things - Groups structure thoughts - Compound representation - SMILES -> Descriptors - SMILES -> Fingerprints # Notes - Meeting for standup at 11am - Remember to use python 3.10 - https://github.com/alan-turing-institute/DSGDec2023Ig - [DeepChem Splitters](https://deepchem.readthedocs.io/en/latest/api_reference/splitters.html) - [OPIG Out-of-distribution generalisation](https://www.blopig.com/blog/2021/06/out-of-distribution-generalisation-and-scaffold-splitting-in-molecular-property-prediction/) ## Group 1 Notes - [Review of Physicochemical Descriptors](https://www.sciencedirect.com/science/article/pii/S0928098720300634) - scaffolding (structurally similar molecules & data leakage) - scaffolding across the three datasets, skikit-learn k-fold - deepchem for scaffolding - three separate models? - sharing scaffolding results & models for comparison between groups 1 & 2 Organisation: - starting to generate PyTorch pipeline - Literature searching - Beginning scaffolding pipeline ## Group 2 Notes - Random forest on the molecular fingerprint as a baseline - Build loader for k-fold cross validation: - First ascertain imbalance of data set - Stratified split for each fold: group split where scaffolds are the groups - Investigate potential for scaffold structuring - DeepChem has library for separating SMILES into scaffolds - Variety of metrics: MCC, F1 score, AUC, Accuracy - Exploration of data sets: bullet point summary - Similarly for existing literature: bullet point summary Data: - Load as NumPy array for larger datasets - Cav: 553 entries, _ (50%) toxic, _ (50%) non-toxic - Herg: 1845 entries, 200 (10%) toxic, 1800 (90%) non-toxic - Note that descriptor data set has 1845 entries, fingerprint data has 1835 - Nav: _ entries, _ (39%) toxic, _ (61%) non-toxic - Data found by running scripts/show-class-inbalance.py - /cav/ignota_turing_dsg_cav_ecfp.csv positives: 278 (0.50%), negatives: 276 (0.50%) - /cav/ignota_turing_dsg_cav_mordred.csv positives: 286 (0.51%), negatives: 276 (0.49%) - /herg/ignota_turing_dsg_herg_ecfp.csv positives: 9124 (0.10%), negatives: 85180 (0.90%) - /herg/ignota_turing_dsg_herg_mordred.csv positives: 9148 (0.10%), negatives: 85279 (0.90%) - /nav/ignota_turing_dsg_nav_ecfp.csv positives: 712 (0.39%), negatives: 1124 (0.61%) - /nav/ignota_turing_dsg_nav_mordred.csv positives: 718 (0.39%), negatives: 1128 (0.61%) - Note that different representations have slightly different counts as some molecules failed to be created and are dropped --- Post Lunch ## Framework Using Pytorch Lightning to setup up pipeline for loading and training data - working on lightning_skeleton branch ## Scaffolding - Use DeepChem ScaffoldSplitter class and / or existing snippet to separate data into scaffolds - Investigate number of scaffolds, size of each scaffold, then see similarity / overlap between scaffolds by UMAP / clustering - Use this to determine: - how train / val / test split should be done - how many buckets we should use, of what size, etc. Nav: - 909 scaffolds, largest of size 99 - 5-fold cross validation, umap plots for each fold - are there disctinct clusters between training and test sets? Or are there overlaps between train-test in chemical space? POA: - fingerprint splitter based on tanimoto similarity - visualise clustering between train and test on UMAP, PCA, & pairwise tanimoto similarity - scaffold splitter - visualise clusering between train and test on UMAP, PCA, & pairwise tanimoto similarity - visualise top n scaffolds in our train and test - maximise distance between tanimoto similarity between scaffolds? Implemented `scaffolding.py` file with `DataSplitter` class. Splits dataset into 5x train-test splits for 5-fold CV. Split data -> train/test & hoholdout Rob: Implemented 'loader.py' that allows loading of 5 fold cross fold based on the 'scaffolding.py' into numpy arrays and into a pytorch Dataset ## Exploratory Data Analysis focused on descriptors rather than fingerprints - correlations between mordred descriptors and toxicity label - notes on mordred descriptors - [list of descriptors](https://mordred-descriptor.github.io/documentation/master/descriptors.html) - [list of constructors and number of descriptors](https://jcheminf.biomedcentral.com/articles/10.1186/s13321-018-0258-y/tables/3) - 2d: 43 constructors - 3d: 5 constructors - 19 descriptors were found to be significant (at 0.01 level) with label - all found significant descriptors are 2D descriptors, constructed by: - ATS (autocorrelation of topological structure) - ATSC (centered ATS) - MATS (Moran coefficient descriptor) - GATS (Gery coeffcient descriptor) - AtomTypeEState (atom type e-state) - number of sCH3 - number of ssNH - number of sOH - sum of sCH3 - sum of ssNH - sum of sOH - ETA delta epsilon (type A) - HBondDonor (hydrogen bond donor descriptor) - RingCount (ring count) - 6-membered aromatic hetero ring count - SLogP (Wildman-Crippen LogP) - [some useful figures to illustrate the structures](https://www.sciencedirect.com/science/article/pii/S0223523420302592) - if time allows: random forest of descriptors - comparison of different representation - [this work](https://pubmed.ncbi.nlm.nih.gov/31138104/) used Dragon (discontinued) for physiochemical descriptors (2432 nonconstant molecular descriptors) of hERG - The neural network model achieved an area under the receiver operating characteristic curve (AUC) of 0.764, with an accuracy of 90.1%, a Matthews correlation coefficient (MCC) of 0.368, a sensitivity of 0.321, and a specificity of 0.967, when ten-fold cross-validation was performed. - its prediction webtool is not available any longer ## Literature Research ### [Images of chemical structures as molecular representations for deep learning](https://link.springer.com/article/10.1557/s43578-022-00628-9) - 2D images of chemical structures Advantages: - allows the user to leverage data augmentation- increases the predictive capabilities of the model by expanding the size and diversity of the datasets. - does not require specialist chemometrics understanding. #### Methodology - Input styles used: - Spectrophore - PubChem Fingerprint - Extended Connectivity Fingerprint (ECFP) - MACCS Key Fingerprint - Mol2Vec - RDKit Descriptors - Mordred Chemical Descriptors (https://github.com/mordred-descriptor/mordred) Rest were generated using - https://github.com/deepchem/deepchem - Data Augmentation: - Solubility data augmentation - Rotations were applied (6 × 30°) to each training example first, and then the resulting images were subjected to three reflections. These reflections occurred along the horizontal axis, along the vertical axis and along both axes together. - Co-crystal data augmentation - Components were additionally stacked in both possible positions, namely component 1 above component 2 and vice versa. Following this, all of the augmentations outlined for the solubility dataset were carried out. ![image_augmentation](https://hackmd.io/_uploads/ryidZDjHT.png) - Model implementation - ResNet architecture was used to evaluate the performance of images as inputs. - Evaluation & Metrics - each dataset was split into training and validation subsets using tenfold cross-validation. The metrics used to assess the performance were accuracy (ACC) and Area Under the Receiver Operating Characteristic Curve (ROC) for the classification task. In the regression task, R2, mean-squared error (MSE), root-mean-squared error (RMSE), and mean absolute error (MAE) were recorded. These metrics were chosen as they are all commonly used in the literature for evaluating AI models. ### Molecular representations in AI-driven drug discovery: a review and practical guide Graph representation - distinction between the notations and file formats built using molecular graphs (e.g. SMILES strings, Molfiles), and the abstract mathematical structure/data structure of a graph itself - general idea: mapping the atoms and bonds that make up a molecule into sets of nodes and edges - Graphs are formally 2D data structures with no spatial relationships between elements (only pair-wise realtionships)-> 3D info must be encoded into a graph representation - nodes: chirality, physiochemical properties, etc - edges: bond lengths, bond types, etc. - all subgraphs are interpretable whereas all SMILES Substrings are not - matrix representations of graphs are node order dependent - graph traversal algorithm determines node order (depth-first or breadth-first search) - Limitations of molecular graph representations - multi-valent bonds (have been introduced by hypergraphs, but not commonly used) - not very compact ### SMILES-BERT - https://dl.acm.org/doi/10.1145/3307339.3342186 ### ChemBERT - pre-trained on ChEMBL database - paper: https://spj.science.org/doi/10.34133/research.0004 - ## Fingerprints (Sarveshwari) Four classes of 2D fingerprint types can be distinguished: (i) dictionary-based, (ii) topological or path-based, (iii) circular fingerprints and (iv) pharmacophores. ### Evaluating fingerprints Just as there is wide variety of fingerprinting algorithms, there are multiple methods for evaluating Virtual Screening performance and little consensus as to which is best. The area under the receiver operating characteristic (**ROC**) curve (**AUC**)/enrichment factor (**EF**) at a given fraction χ of the data set * Advantage: The advantage of AUC is that it is bounded, running from 0 to 1 with 0.5 corresponding to randomness, and that it is independent of the ratio of actives to inactives and other external parameters. * Disadvantage: The AUC method has been critized as being inappropriate for comparing VS methods as it is not sufficiently sensitive to early recognition. The EF explicitly measures early recognition but it is dependent on the ratio of actives to inactives and the choice of χ. To try and overcome these limitations numerous other evaluation methods, such as robust initial enhancement (**RIE**) and Boltzmann-enhanced discrimination of ROC (**BEDROC**), have been proposed. * Advantage: The RIE uses a continuously decreasing exponential weight as a function of rank and is thus sensitive to early recognition. * Disadvantage: It is dependent both on an adjustable parameter, the exponential weight, and the ratio of actives to inactives. RIE values can therefore not be easily compared between different data sets. The BEDROC is constructed by, in essence, forcing the RIE to be bounded by 0 and 1, avoiding the dependence on the active/inactive ratio. | Name | Type | |---------- | ----------| | ECFP4 | circular fingerprints with bit string | | | | long ECFP4 | circular fingerprints with bit string | | | ECFP6 | circular fingerprints with bit string | | | | long ECFP6 | circular fingerprints with bit string | | | | FCFP4 | circular fingerprints with bit string | | | | FCFC4 | circular fingerprints with count vector | | | | ECFC4 | circular fingerprints with count vector | | | | MACCS | dictionary-based fingerprint | | | | RDK5 | Path-based fingerprint | | | | long Avalon | Path-based fingerprint | | | | Avalon | Path-based fingerprint | | | | TT | Path-based fingerprint | | | | AP | Path-based fingerprint | | | | ECFC0 | Text | Text | | #### ECFPs ECFP belongs to the best performing fingerprints in small molecule virtual screening and target prediction benchmarks, together with the related MinHashed fingerprint MHFP6. ECFPs are circular fingerprints with a number of useful qualities: they can be very **rapidly calculated**; they are **not predefined** and can represent an essentially infinite number of different molecular features (including stereochemical information); their **features represent the presence of particular substructures**, allowing easier interpretation of analysis results; and the ECFP algorithm can be **tailored to generate different types of circular fingerprints**, optimized for different uses. *Being already formatted as a vector with numerical entries, they are well suited to be used as features in ML models for chemistry.* **Comparison to other fingerprint methods** ECFPs have several advantages when contrasted to ***predefined substructural keys***. 1. Since the set of keys is fixed, the system may not contain keys appropriate to the novel structural variation in a given compound library. 2. The keys are designed for substructure search, which limits their utility for activity modeling and categorization (see the Hydrogen-Filled Substructural Features Section for this discussion). 3. Even the most intelligently selected set of keys18 is limited in coverage, making it unlikely that more detailed and specific substructures particular to some unusual activity class would be represented. Another commonly used class of fingerprints is available through Daylight Chemical Information Systems. It uses features based upon the presence of paths of varying lengths containing specific atom types. This generates a sparse binary vector, which is commonly “folded” to a bitset of specified size (e.g., 1024 bits) to reduce its size for ease of manipulation. ECFPs have several advantages when contrasted to such ***pathbased fingerprints***. 1. Daylight keys were designed for substructure and similarity search, which limits their effectiveness for activity modeling. If some path feature is found to be useful during an analysis, then interpretation can still be difficult, as a given path may not correspond to an easily identifiable function group or substructural entity. 2. Further complicated but often relevant chemical attachment patterns, such as quaternary centers or cyclic patterns, cannot be directly represented by a path-based description. Path-based schemes also cannot capture atom-based stereochemistry. **Disadvantages of ECFPs (ECFP4) and MHFP6:** Both have a poor perception of the global features of molecules such as size and shape. They also fail at perceiving structural differences that may be important in larger molecules, such as distinguishing between regioisomers in extended ring systems (e.g. 2,7- versus 2,8-dichlorodioxin), between linkers of different lengths, or between scrambled peptide sequences of identical composition and length. ### Contemporary fingerprints #### MAP4 (MinHashed Atom-Pair fingerprint up to four bonds) The ECFP4 and MHFP6 limitations can be addressed by using atom-pair fingerprints, which encode molecular shape and are often used for scaffold-hopping. However, they do not encode molecular structure in detail and perform poorly in small molecule benchmarking studies compared to substructure fingerprints such as ECFP4 and MHFP6. MAP4 encodes atom pairs and their bond distance similarly to the **AP fingerprint implemented by RDKit**, however in MAP4 atom characteristics are replaced by the circular substructure around each atom of the pair, written in **SMILES format**. MAP4 uses the same MinHashing technique as **MHFP6**, a principle borrowed from natural language processing which enables fast similarity searches in very large databases by locality sensitive hashing (LSH). LSH is a technique that allows the creation of self-tuning indexes, which are then used to generate a forest of trees that can be traversed for an approximate but fast similarity search. --- ## Alternative representations (Sarveshwari) * **Coulomb Matrix** Ab-initio electronic structure calculations typically require a set of nuclear charges {Z} and the corresponding Cartesian coordinates {R} as input. The Coulomb Matrix (CM) M, proposed by Rupp et al., encodes this information by use of the atomic self-energies and internuclear Coulomb repulsion operator. In the construction of coulomb matrix, we first use the nuclear charges and distance matrix generated by RDKit to acquire the original coulomb matrix, then an optional random atom index sorting and binary expansion transformation can be applied during training in order to achieve atom index invariance, as reported by Montavon et al. *Example code for converting from SMILES to Coulomb Matrices **using RDKit in Python**:* ``` import rdkit from rdkit import Chem from rdkit.Chem import AllChem from rdkit.Chem.rdMolDescriptors import CalcCoulombMat smiles = 'CN1C=NC2=C1C(=O)N(C(=O)N2C)C' caffeine = Chem.MolFromSmiles(smiles, sanitize=True) caffeine = Chem.AddHs(caffeine) AllChem.EmbedMolecule(caffeine) AllChem.UFFOptimizeMolecule(caffeine) CM = CalcCoulombMat(caffeine) list(CM[0]) ``` *The implementation of coulomb matrix descriptors by Matthias Rupp et. al. **using the ChemML library*** ``` from chemml.chem import CoulombMatrix, Molecule m1 = Molecule('c1ccc1', 'smiles') m2 = Molecule('CNC', 'smiles') m3 = Molecule('CC', 'smiles') m4 = Molecule('CCC', 'smiles') molecules = [m1, m2, m3, m4] for mol in molecules: mol.to_xyz(optimizer='UFF') cm = CoulombMatrix(cm_type='SC', n_jobs=-1) features = cm.represent(molecules) ``` In the above, cm_type is Coulomb matrix type, which can be one of the five below - * ‘Unsorted_Matrix’ or ‘UM’ * ‘Unsorted_Triangular’ or ‘UT’ * ‘Eigenspectrum’ or ‘E’ * ‘Sorted_Coulomb’ or ‘SC’ (default) * ‘Random_Coulomb’ or ‘RC’ *Generate Coulomb matrices **from scratch using RDKit**:* ``` """ Generate coulomb matrices for molecules. See Montavon et al., _New Journal of Physics_ __15__ (2013) 095003. """ import numpy as np from typing import Any, List, Optional from deepchem.utils.typing import RDKitMol from deepchem.utils.data_utils import pad_array from deepchem.feat.base_classes import MolecularFeaturizer class CoulombMatrix(MolecularFeaturizer): """Calculate Coulomb matrices for molecules. Coulomb matrices provide a representation of the electronic structure of a molecule. For a molecule with `N` atoms, the Coulomb matrix is a `N X N` matrix where each element gives the strength of the electrostatic interaction between two atoms. The method is described in more detail in [1]_. Examples -------- >>> import deepchem as dc >>> featurizers = dc.feat.CoulombMatrix(max_atoms=23) >>> input_file = 'deepchem/feat/tests/data/water.sdf' # really backed by water.sdf.csv >>> tasks = ["atomization_energy"] >>> loader = dc.data.SDFLoader(tasks, featurizer=featurizers) >>> dataset = loader.create_dataset(input_file) References ---------- .. [1] Montavon, Grégoire, et al. "Learning invariant representations of molecules for atomization energy prediction." Advances in neural information processing systems. 2012. Note ---- This class requires RDKit to be installed. """ def __init__(self, max_atoms: int, remove_hydrogens: bool = False, randomize: bool = False, upper_tri: bool = False, n_samples: int = 1, seed: Optional[int] = None): """Initialize this featurizer. Parameters ---------- max_atoms: int The maximum number of atoms expected for molecules this featurizer will process. remove_hydrogens: bool, optional (default False) If True, remove hydrogens before processing them. randomize: bool, optional (default False) If True, use method `randomize_coulomb_matrices` to randomize Coulomb matrices. upper_tri: bool, optional (default False) Generate only upper triangle part of Coulomb matrices. n_samples: int, optional (default 1) If `randomize` is set to True, the number of random samples to draw. seed: int, optional (default None) Random seed to use. """ self.max_atoms = int(max_atoms) self.remove_hydrogens = remove_hydrogens self.randomize = randomize self.upper_tri = upper_tri self.n_samples = n_samples if seed is not None: seed = int(seed) self.seed = seed def _featurize(self, datapoint: RDKitMol, **kwargs) -> np.ndarray: """ Calculate Coulomb matrices for molecules. If extra randomized matrices are generated, they are treated as if they are features for additional conformers. Since Coulomb matrices are symmetric, only the (flattened) upper triangular portion is returned. Parameters ---------- datapoint: rdkit.Chem.rdchem.Mol RDKit Mol object Returns ------- np.ndarray The coulomb matrices of the given molecule. The default shape is `(num_confs, max_atoms, max_atoms)`. If num_confs == 1, the shape is `(max_atoms, max_atoms)`. """ if 'mol' in kwargs: datapoint = kwargs.get("mol") raise DeprecationWarning( 'Mol is being phased out as a parameter, please pass "datapoint" instead.' ) features = self.coulomb_matrix(datapoint) if self.upper_tri: features = np.asarray( [f[np.triu_indices_from(f)] for f in features]) if features.shape[0] == 1: # `(1, max_atoms, max_atoms)` -> `(max_atoms, max_atoms)` features = np.squeeze(features, axis=0) return features def coulomb_matrix(self, mol: RDKitMol) -> np.ndarray: """ Generate Coulomb matrices for each conformer of the given molecule. Parameters ---------- mol: rdkit.Chem.rdchem.Mol RDKit Mol object Returns ------- np.ndarray The coulomb matrices of the given molecule """ try: from rdkit import Chem from rdkit.Chem import AllChem except ModuleNotFoundError: raise ImportError("This class requires RDKit to be installed.") # Check whether num_confs >=1 or not num_confs = len(mol.GetConformers()) if num_confs == 0: mol = Chem.AddHs(mol) AllChem.EmbedMolecule(mol, AllChem.ETKDG()) if self.remove_hydrogens: mol = Chem.RemoveHs(mol) n_atoms = mol.GetNumAtoms() z = [atom.GetAtomicNum() for atom in mol.GetAtoms()] rval = [] for conf in mol.GetConformers(): d = self.get_interatomic_distances(conf) m = np.outer(z, z) / d m[range(n_atoms), range(n_atoms)] = 0.5 * np.array(z)**2.4 if self.randomize: for random_m in self.randomize_coulomb_matrix(m): random_m = pad_array(random_m, self.max_atoms) rval.append(random_m) else: m = pad_array(m, self.max_atoms) rval.append(m) return np.asarray(rval) def randomize_coulomb_matrix(self, m: np.ndarray) -> List[np.ndarray]: """Randomize a Coulomb matrix as decribed in [1]_: 1. Compute row norms for M in a vector row_norms. 2. Sample a zero-mean unit-variance noise vector e with dimension equal to row_norms. 3. Permute the rows and columns of M with the permutation that sorts row_norms + e. Parameters ---------- m: np.ndarray Coulomb matrix. Returns ------- List[np.ndarray] List of the random coulomb matrix References ---------- .. [1] Montavon et al., New Journal of Physics, 15, (2013), 095003 """ rval = [] row_norms = np.asarray([np.linalg.norm(row) for row in m], dtype=float) rng = np.random.RandomState(self.seed) for i in range(self.n_samples): e = rng.normal(size=row_norms.size) p = np.argsort(row_norms + e) new = m[p][:, p] # permute rows first, then columns rval.append(new) return rval @staticmethod def get_interatomic_distances(conf: Any) -> np.ndarray: """ Get interatomic distances for atoms in a molecular conformer. Parameters ---------- conf: rdkit.Chem.rdchem.Conformer Molecule conformer. Returns ------- np.ndarray The distances matrix for all atoms in a molecule """ n_atoms = conf.GetNumAtoms() coords = [ # Convert AtomPositions from Angstrom to bohr (atomic units) conf.GetAtomPosition(i).__idiv__(0.52917721092) for i in range(n_atoms) ] d = np.zeros((n_atoms, n_atoms), dtype=float) for i in range(n_atoms): for j in range(i): d[i, j] = coords[i].Distance(coords[j]) d[j, i] = d[i, j] return d class CoulombMatrixEig(CoulombMatrix): """Calculate the eigenvalues of Coulomb matrices for molecules. This featurizer computes the eigenvalues of the Coulomb matrices for provided molecules. Coulomb matrices are described in [1]_. Examples -------- >>> import deepchem as dc >>> featurizers = dc.feat.CoulombMatrixEig(max_atoms=23) >>> input_file = 'deepchem/feat/tests/data/water.sdf' # really backed by water.sdf.csv >>> tasks = ["atomization_energy"] >>> loader = dc.data.SDFLoader(tasks, featurizer=featurizers) >>> dataset = loader.create_dataset(input_file) References ---------- .. [1] Montavon, Grégoire, et al. "Learning invariant representations of molecules for atomization energy prediction." Advances in neural information processing systems. 2012. """ def __init__(self, max_atoms: int, remove_hydrogens: bool = False, randomize: bool = False, n_samples: int = 1, seed: Optional[int] = None): """Initialize this featurizer. Parameters ---------- max_atoms: int The maximum number of atoms expected for molecules this featurizer will process. remove_hydrogens: bool, optional (default False) If True, remove hydrogens before processing them. randomize: bool, optional (default False) If True, use method `randomize_coulomb_matrices` to randomize Coulomb matrices. n_samples: int, optional (default 1) If `randomize` is set to True, the number of random samples to draw. seed: int, optional (default None) Random seed to use. """ self.max_atoms = int(max_atoms) self.remove_hydrogens = remove_hydrogens self.randomize = randomize self.n_samples = n_samples if seed is not None: seed = int(seed) self.seed = seed def _featurize(self, datapoint: RDKitMol, **kwargs) -> np.ndarray: """ Calculate eigenvalues of Coulomb matrix for molecules. Eigenvalues are returned sorted by absolute value in descending order and padded by max_atoms. Parameters ---------- datapoint: rdkit.Chem.rdchem.Mol RDKit Mol object Returns ------- np.ndarray The eigenvalues of Coulomb matrix for molecules. The default shape is `(num_confs, max_atoms)`. If num_confs == 1, the shape is `(max_atoms,)`. """ if 'mol' in kwargs: datapoint = kwargs.get("mol") raise DeprecationWarning( 'Mol is being phased out as a parameter, please pass "datapoint" instead.' ) cmat = self.coulomb_matrix(datapoint) features_list = [] for f in cmat: w, v = np.linalg.eig(f) w_abs = np.abs(w) sortidx = np.argsort(w_abs) sortidx = sortidx[::-1] w = w[sortidx] f = pad_array(w, self.max_atoms) features_list.append(f) features = np.asarray(features_list) if features.shape[0] == 1: # `(1, max_atoms)` -> `(max_atoms,)` features = np.squeeze(features, axis=0) return features ``` * **Grid Featurizer** The grid featurizer was inspired by the NNscore featurizer and SPLIF69 but optimized for speed, robustness, and generalizability. The intermolecular interactions enumerated by the featurizer include salt bridges and hydrogen bonding between protein and ligand, intraligand circular fingerprints, intra-protein circular fingerprints, and protein-ligand SPLIF fingerprints. RdkitGridFeaturizer featurizes protein-ligand complex using flat features or a 3D grid (in which each voxel is described with a vector of features). ``` ``` * **Symmetry Function** Symmetry function, first introduced by Behler and Parrinello, is another common encoding of atomic coordinates information. It focuses on preserving the rotational and permutation symmetry of the system. The local environment of an atom in the molecule is expressed as a series of radial and angular symmetry functions with different distance and angle cutoffs, the former focusing on distances between atom pairs and the latter focusing on angles formed within triplets of atoms. In the paper titled "ANI-1: An extensible neural network potential with DFT accuracy at force field computational cost", the authors demonstrate how a deep neural network (NN) trained on quantum mechanical (QM) DFT calculations can learn an accurate and fully transferable potential for organic molecules. Authors introduce ANAKIN-ME (Accurate NeurAl networK engINe for Molecular Energies) or ANI in short, which can be implemented using the DeepChem package in Python using the following command: ``` deepchem.models.ANIRegression ``` Dataset preparation from SMILES to symmetry function in Python: ``` import numpy as np from deepchem.utils.typing import RDKitMol from deepchem.utils.data_utils import pad_array from deepchem.feat.base_classes import MolecularFeaturizer from deepchem.feat.molecule_featurizers.atomic_coordinates import AtomicCoordinates class BPSymmetryFunctionInput(MolecularFeaturizer): """Calculate symmetry function for each atom in the molecules This method is described in [1]_. Examples -------- >>> import deepchem as dc >>> smiles = ['C1C=CC=CC=1'] >>> featurizer = dc.feat.BPSymmetryFunctionInput(max_atoms=10) >>> features = featurizer.featurize(smiles) >>> type(features[0]) <class 'numpy.ndarray'> >>> features[0].shape # (max_atoms, 4) (10, 4) References ---------- .. [1] Behler, Jörg, and Michele Parrinello. "Generalized neural-network representation of high-dimensional potential-energy surfaces." Physical review letters 98.14 (2007): 146401. Note ---- This class requires RDKit to be installed. """ def __init__(self, max_atoms: int): """Initialize this featurizer. Parameters ---------- max_atoms: int The maximum number of atoms expected for molecules this featurizer will process. """ self.max_atoms = max_atoms self.coordfeat = AtomicCoordinates(use_bohr=True) def _featurize(self, datapoint: RDKitMol, **kwargs) -> np.ndarray: """Calculate symmetry function. Parameters ---------- datapoint: rdkit.Chem.rdchem.Mol RDKit Mol object Returns ------- np.ndarray A numpy array of symmetry function. The shape is `(max_atoms, 4)`. """ if 'mol' in kwargs: datapoint = kwargs.get("mol") raise DeprecationWarning( 'Mol is being phased out as a parameter, please pass "datapoint" instead.' ) coordinates = self.coordfeat._featurize(datapoint) atom_numbers = np.array( [atom.GetAtomicNum() for atom in datapoint.GetAtoms()]) atom_numbers = np.expand_dims(atom_numbers, axis=1) assert atom_numbers.shape[0] == coordinates.shape[0] features = np.concatenate([atom_numbers, coordinates], axis=1) return pad_array(features, (self.max_atoms, 4)) ``` * **Weave** Weave featurization encodes both local chemical environment and connectivity of atoms in a molecule. Atomic feature vectors are exactly the same, while connectivity is represented by more detailed pair features instead of neighbor listing. The weave featurization calculates a feature vector for each pair of atoms in the molecule, including bond properties (if directly connected), graph distance and ring info, forming a feature matrix. The method supports graph based models that utilize properties of both nodes (atoms) and edges (bonds). Weave convolutions were introduced in Kearnes, Steven, et al. Unlike Duvenaud graph convolutions, weave convolutions require a quadratic matrix of interaction descriptors for each pair of atoms. These extra descriptors may provide for additional descriptive power but at the cost of a larger featurized dataset. ``` import deepchem as dc mols = ["CCC"] featurizer = dc.feat.WeaveFeaturizer() features = featurizer.featurize(mols) type(features[0]) features[0].get_num_atoms() # 3 atoms in compound features[0].get_num_features() # feature size type(features[0].get_atom_features()) features[0].get_atom_features().shape type(features[0].get_pair_features()) features[0].get_pair_features().shape ``` * **Image** *SMILEs to image code:* ``` """ Featurizer implementations used in ChemCeption models. SmilesToImage featurizer for ChemCeption models taken from https://arxiv.org/abs/1710.02238 """ import numpy as np from deepchem.utils.typing import RDKitMol from deepchem.feat.base_classes import MolecularFeaturizer class SmilesToImage(MolecularFeaturizer): """Convert SMILES string to an image. SmilesToImage Featurizer takes a SMILES string, and turns it into an image. Details taken from [1]_. The default size of for the image is 80 x 80. Two image modes are currently supported - std & engd. std is the gray scale specification, with atomic numbers as pixel values for atom positions and a constant value of 2 for bond positions. engd is a 4-channel specification, which uses atom properties like hybridization, valency, charges in addition to atomic number. Bond type is also used for the bonds. The coordinates of all atoms are computed, and lines are drawn between atoms to indicate bonds. For the respective channels, the atom and bond positions are set to the property values as mentioned in the paper. References ---------- .. [1] Goh, Garrett B., et al. "Using rule-based labels for weak supervised learning: a ChemNet for transferable chemical property prediction." Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining. 2018. Notes ----- This class requires RDKit to be installed. """ def __init__(self, img_size: int = 80, res: float = 0.5, max_len: int = 250, img_spec: str = "std"): """ Parameters ---------- img_size: int, default 80 Size of the image tensor res: float, default 0.5 Displays the resolution of each pixel in Angstrom max_len: int, default 250 Maximum allowed length of SMILES string img_spec: str, default std Indicates the channel organization of the image tensor """ if img_spec not in ["std", "engd"]: raise ValueError( "Image mode must be one of std or engd. {} is not supported".format( img_spec)) self.img_size = img_size self.max_len = max_len self.res = res self.img_spec = img_spec self.embed = int(img_size * res / 2) def _featurize(self, mol: RDKitMol) -> np.ndarray: """Featurizes a single SMILE into an image. Parameters ---------- mol: rdkit.Chem.rdchem.Mol RDKit Mol object Returns ------- np.ndarray A 3D array of image, the shape is `(img_size, img_size, 1)`. If the length of SMILES is longer than `max_len`, this value is an empty array. """ try: from rdkit import Chem from rdkit.Chem import AllChem except ModuleNotFoundError: raise ImportError("This class requires RDKit to be installed.") smile = Chem.MolToSmiles(mol) if len(smile) > self.max_len: return np.array([]) cmol = Chem.Mol(mol.ToBinary()) cmol.ComputeGasteigerCharges() AllChem.Compute2DCoords(cmol) atom_coords = cmol.GetConformer(0).GetPositions() if self.img_spec == "std": # Setup image img = np.zeros((self.img_size, self.img_size, 1)) # Compute bond properties bond_props = np.array( [[2.0, bond.GetBeginAtomIdx(), bond.GetEndAtomIdx()] for bond in mol.GetBonds()]) # Compute atom properties atom_props = np.array([[atom.GetAtomicNum()] for atom in cmol.GetAtoms()]) bond_props = bond_props.astype(np.float32) atom_props = atom_props.astype(np.float32) else: # Setup image img = np.zeros((self.img_size, self.img_size, 4)) # Compute bond properties bond_props = np.array([[ bond.GetBondTypeAsDouble(), bond.GetBeginAtomIdx(), bond.GetEndAtomIdx() ] for bond in mol.GetBonds()]) # Compute atom properties atom_props = np.array([[ atom.GetAtomicNum(), atom.GetProp("_GasteigerCharge"), atom.GetExplicitValence(), atom.GetHybridization().real, ] for atom in cmol.GetAtoms()]) bond_props = bond_props.astype(np.float32) atom_props = atom_props.astype(np.float32) partial_charges = atom_props[:, 1] if np.any(np.isnan(partial_charges)): return np.array([]) frac = np.linspace(0, 1, int(1 / self.res * 2)) # Reshape done for proper broadcast frac = frac.reshape(-1, 1, 1) bond_begin_idxs = bond_props[:, 1].astype(int) bond_end_idxs = bond_props[:, 2].astype(int) # Reshapes, and axes manipulations to facilitate vector processing. begin_coords = atom_coords[bond_begin_idxs] begin_coords = np.expand_dims(begin_coords.T, axis=0) end_coords = atom_coords[bond_end_idxs] end_coords = np.expand_dims(end_coords.T, axis=0) # Draw a line between the two atoms. # The coordinates of this line, are indicated in line_coords line_coords = frac * begin_coords + (1 - frac) * end_coords # Turn the line coordinates into image positions bond_line_idxs = np.ceil( (line_coords[:, 0] + self.embed) / self.res).astype(int) bond_line_idys = np.ceil( (line_coords[:, 1] + self.embed) / self.res).astype(int) # Set the bond line coordinates to the bond property used. img[bond_line_idxs, bond_line_idys, 0] = bond_props[:, 0] # Turn atomic coordinates into image positions atom_idxs = np.round( (atom_coords[:, 0] + self.embed) / self.res).astype(int) atom_idys = np.round( (atom_coords[:, 1] + self.embed) / self.res).astype(int) # Set the atom positions in image to different atomic properties in channels img[atom_idxs, atom_idys, :] = atom_props return img ``` --- ### Scaffold Splitting - - dgl-lifesci - https://github.com/awslabs/dgl-lifesci - GNN-based modeling on custom datasets for molecular property prediction, reaction prediction, and molecule generation. - ### A systematic study of key elements underlying molecular property prediction - MoleculeNet ## References *Fingerprint Representations* * Rogers, D., & Hahn, M. (2010). Extended-Connectivity Fingerprints. Journal of Chemical Information and Modeling, 50(5), 742–754. doi:10.1021/ci100050t * https://docs.chemaxon.com/display/docs/extended-connectivity-fingerprint-ecfp.md * Capecchi, A., Probst, D. & Reymond, JL. One molecular fingerprint to rule them all: drugs, biomolecules, and the metabolome. J Cheminform 12, 43 (2020). * Riniker, S., Landrum, G.A. Open-source platform to benchmark fingerprints for ligand-based virtual screening. J Cheminform 5, 26 (2013). https://doi.org/10.1186/1758-2946-5-26 * Wu, Zhenqin et al. “MoleculeNet: a benchmark for molecular machine learning.” Chemical science vol. 9,2 513-530. 31 Oct. 2017 * Goh, Garrett B., et al. “Using rule-based labels for weak supervised learning: a ChemNet for transferable chemical property prediction.” Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining. 2018. * Kearnes, Steven, et al. “Molecular graph convolutions: moving beyond fingerprints.” Journal of computer-aided molecular design 30.8 (2016): 595-608. *Graph Molecular Representations* 1. Xiong, Z. et al. Pushing the boundaries of molecular representation for drug discovery with the graph attention mechanism. J. Med. Chem. 63, 8749–8760 (2019) 2. Na, G. S., Chang, H. & Kim, H. W. Machine-guided representation for accurate graph-based molecular machine learning. Phys. Chem. Chem. Phys. 22, 18526–18535 (2020) 3. Jiang, D. et al. Could graph neural networks learn better molecular representation for drug discovery? A comparison study of descriptor-based and graph-based models. J. Cheminformatics 13, 1–23 (2021) 4. Wilkinson, M.R., Martinez-Hernandez, U., Wilson, C.C. et al. Images of chemical structures as molecular representations for deep learning. Journal of Materials Research 37, 2293–2303 (2022). https://doi.org/10.1557/s43578-022-00628-9 *Molecular property prediction* 1. Deng, J., Yang, Z., Wang, H. et al. A systematic study of key elements underlying molecular property prediction. Nat Commun 14, 6395 (2023). https://doi.org/10.1038/s41467-023-41948-6 ## Baseline models - Based on tabular data - dropped smiles, nans - Models: - Random forest **Random split stratified by class ** ECFP | | CAV | HERG | NAV | | --------- | ---- | ---- | ---- | | Accuracy: |0.67(+-0.11)|0.91 (+-0.00) |0.6(+-0.06)| |F1: |0.64(+-0.13)| 0.88 (+-0.00) |0.55(+-0.04)| | MCC: | 0.39(+-0.22) | 0.29 (+-0.01) | 0.11(+-0.12) | MORDRED | | CAV | HERG | NAV | | --- | ---- | ---- | ---- | | Accuracy: | Text | Text | Text | Scaffold split: ECFP | | CAV | HERG | NAV | | --------- | ---- | ---- | ---- | | Accuracy: |0.72(+-0.04)|0.91 (+-0.00) |0.74(+-0.04)| |F1: |0.71(+-0.06)|0.87 (+-0.00) |0.71(+-0.06)| | MCC: | 0.44(+-0.11) | 0.25 (+-0.03) | 0.44(+-0.11) | MORDRED | | CAV | HERG | NAV | | --------- | ---- | ---- | ---- | | Accuracy:| 0.77 (+-0.04) | | 0.77 (+-0.04) | |F1: | 0.77 (+-0.05) | | 0.77 (+-0.05) |MCC: | 0.52 (+-0.10) | Text | 0.52 (+-0.10) ADABOOST [CM sum across folds]: * NAV, ECFP: | Predicted: | 0 | 1 | | -------- | -------- | -------- | |0|942| 181| | 1 | 290| 422| * NAV, MORDRED: | Predicted: | 0 | 1 | | -------- | -------- | -------- | |0|938 |189| | 1 | 304 |414 | * CAV, ECFP: | Predicted: | 0 | 1 | | -------- | -------- | -------- | |0|218 | 57| | 1 | 75 |203| * CAV, MORDRED: | Predicted: | 0 | 1 | | -------- | -------- | -------- | |0|215 | 60| | 1 | 59| 227| * HERG, ECFP: | Predicted: | 0 | 1 | | -------- | -------- | -------- | |0|83548 |1631 | 1 | 5874 | 3250| * HERG, MORDRED: | Predicted: | 0 | 1 | | -------- | -------- | -------- | |0|83372 | 1906| | 1 | 4781 | 4367|

Syntax	Example	Reference
# Header	Header	基本排版
- Unordered List	Unordered List
1. Ordered List	Ordered List
- [ ] Todo List	Todo List
> Blockquote	Blockquote
Bold font	Bold font
Italics font	Italics font
~~Strikethrough~~	~~Strikethrough~~
19^th^	19^th
H~2~O	H₂O
++Inserted text++	Inserted text
==Marked text==	Marked text
[link text](https:// "title")	Link
![image alt](https:// "title")	Image
`Code`	`Code`	在筆記中貼入程式碼
```javascript var i = 0; ```	`var i = 0;`
:smile:		Emoji list
{%youtube youtube_id %}	Externals
$L^aT_eX$	L^aT_eX
:::info This is a alert area. :::	This is a alert area.