Making features for the building blocks SMILES (Physicochemical Descriptors, Fingerprints)\n
Similarity comparisons and visualisation (multi-dimensional scaling, heatmaps)\n
\n The same approach can be used to generate the features for the train and test molecules (molecule_smiles). \n In this notebook I use and cite the ChemmineR: Cheminformatics Toolkit for R tutorial and add-on package fmcsR\n
","metadata":{}},{"cell_type":"code","source":"suppressMessages({ \n if (!requireNamespace(\"BiocManager\", quietly=TRUE))\n install.packages(\"BiocManager\")\n BiocManager::install(\"ChemmineR\", update = FALSE, quiet = TRUE)\n BiocManager::install(\"fmcsR\", update = FALSE, quiet = TRUE) \n library(ChemmineR)\n library(fmcsR)\n library(data.table)\n library(qs)\n library(foreach)\n library(pheatmap)\n library(ggplot2)\n})","metadata":{"_uuid":"051d70d956493feee0c6d64651c6a088724dca2a","_execution_state":"idle","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Helper functions**","metadata":{}},{"cell_type":"code","source":"shape <- function(dt) {\n cat(\"\\nShape: \", ncol(dt), \" columns x \",\n nrow(dt), \" rows\", sep = \"\")\n}\n\nfig <- function(width, heigth) {\n options(repr.plot.width = width, repr.plot.height = heigth)\n}","metadata":{"execution":{"iopub.status.busy":"2024-04-09T15:26:30.793541Z","iopub.execute_input":"2024-04-09T15:26:30.823359Z","iopub.status.idle":"2024-04-09T15:26:30.838475Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"****","metadata":{}},{"cell_type":"markdown","source":"# Load the data","metadata":{}},{"cell_type":"markdown","source":"Here I load all the unique SMILES from all three blocks (BB) from the test and train sets (extracted in [this notebook](https://www.kaggle.com/code/antoninadolgorukova/belka-reading-and-quick-stats)).","metadata":{}},{"cell_type":"code","source":"bb_smiles <- fread(\"/kaggle/input/belka-supplementary-calcs-and-data-for-ml/smiles/all_bb_smiles_by_bb.csv\")\nbb_smiles[, .N, by = c(\"set\", \"BB\")]","metadata":{"execution":{"iopub.status.busy":"2024-04-09T15:26:30.841125Z","iopub.execute_input":"2024-04-09T15:26:30.842483Z","iopub.status.idle":"2024-04-09T15:26:30.913489Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"I also create an aggregated table by building block (BB) and data set (train, test, or both) to use when plotting.","metadata":{}},{"cell_type":"code","source":"bb_smiles_agg <- bb_smiles[\n , .(BB = paste(unique(BB), collapse = \", \"),\n set = paste(unique(set), collapse = \"/\")),\n by = \"smiles\"]\n\nbb_smiles_agg[, .N, by = c(\"BB\", \"set\")][order(BB)]","metadata":{"execution":{"iopub.status.busy":"2024-04-09T16:02:43.997428Z","iopub.execute_input":"2024-04-09T16:02:44.003530Z","iopub.status.idle":"2024-04-09T16:02:44.090927Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n 💡 Note: The table summarizes the counts of unique smiles for each building block (BB) within different datasets (set). \n \n
\n","metadata":{}},{"cell_type":"markdown","source":"The SMILES of all blocks were converted to sfd format on my PC and stored in the [BELKA: supplementary calcs & data for ML dataset](https://www.kaggle.com/datasets/antoninadolgorukova/belka-supplementary-calcs-and-data-for-ml/data), because the required package [ChemmineOB](https://bioconductor.org/packages/release/bioc/html/ChemmineOB.html) and the OpenBabel software are not available in kagge. The code to convert smiles to sdf is in the commented section below. ","metadata":{}},{"cell_type":"markdown","source":"
\n 💡 Note: \n
\n
In chemistry and cheminformatics, SDF stands for Structure Data File. An SDF file is a standard file format used to represent chemical structures and associated data. It's a text-based format that stores information about molecules, including atoms, bonds, coordinates, and various properties.\n
The molecules have IDs like 'CMP'+ number, and are in the same order as in the characted vector of SMILES strings bb_smiles.\n
\n
\n\n","metadata":{}},{"cell_type":"code","source":"bb_smiles_sdfset <- read.SDFset(\"/kaggle/input/belka-supplementary-calcs-and-data-for-ml/smiles/bb_smiles_sdfset.sdf\")\nbb_smiles_sdfset\n\ncat(\"\\nSingle molecule from SDFset:\\n\")\nas(bb_smiles_sdfset[[1]], \"list\")","metadata":{"execution":{"iopub.status.busy":"2024-04-09T15:26:30.915970Z","iopub.execute_input":"2024-04-09T15:26:30.917279Z","iopub.status.idle":"2024-04-09T15:26:33.720942Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"You can visualize compound structures with a standard web browser using the [ChemMine Tools](https://chemminetools.ucr.edu/) online service. This line ran on PC generate link and open in ChemMine Tools in your browser.","metadata":{}},{"cell_type":"code","source":"# sdf.visualize(bb_smiles_sdfset[1:5])","metadata":{"execution":{"iopub.status.busy":"2024-04-09T15:26:33.723788Z","iopub.execute_input":"2024-04-09T15:26:33.725274Z","iopub.status.idle":"2024-04-09T15:26:33.738096Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Make a dictionary with smiles and ids in SDF**","metadata":{}},{"cell_type":"code","source":"smiles_dict <- data.table(\n smiles = sapply(1:length(bb_smiles_sdfset), function(i) {\n bb_smiles_sdfset[[i]]@header[\"Molecule_Name\"]\n }),\n id = cid(bb_smiles_sdfset)\n)\nhead(smiles_dict)","metadata":{"execution":{"iopub.status.busy":"2024-04-09T15:26:33.742187Z","iopub.execute_input":"2024-04-09T15:26:33.743803Z","iopub.status.idle":"2024-04-09T15:26:33.805932Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Generate Physicochemical Descriptors","metadata":{}},{"cell_type":"markdown","source":"Basic compound descriptors such as molecular formula (MF), molecular weight (MW), and atomic and functional group frequencies can be calculated using ChemmineR functions.","metadata":{}},{"cell_type":"code","source":"mol_features <- data.table(\n SMILES = smiles_dict$smiles,\n MF = MF(bb_smiles_sdfset, addH = TRUE),\n MW = MW(bb_smiles_sdfset, addH = TRUE),\n Ncharges = sapply(bonds(bb_smiles_sdfset, type = \"charge\"), length),\n atomcountMA(bb_smiles_sdfset, addH = TRUE),\n groups(bb_smiles_sdfset, type = \"countMA\"),\n rings(bb_smiles_sdfset, type = \"count\", arom = \"TRUE\"))\n\nhead(mol_features)\nshape(mol_features)","metadata":{"execution":{"iopub.status.busy":"2024-04-09T15:26:33.809978Z","iopub.execute_input":"2024-04-09T15:26:33.811397Z","iopub.status.idle":"2024-04-09T15:27:04.626686Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"names(mol_features)","metadata":{"execution":{"iopub.status.busy":"2024-04-09T15:27:04.630531Z","iopub.execute_input":"2024-04-09T15:27:04.631965Z","iopub.status.idle":"2024-04-09T15:27:04.650042Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"About rings: The function `rings` identifies all possible rings in molecules using the exhaustive ring perception algorithm from Hanser et al (1996). In addition, the function can return all smallest possible rings as well as aromaticity information for each ring. A ring is considered aromatic if it meets the following requirements: \n- all atoms in the ring need to be sp2 hybridized. This means each atom has to have a double bond or at least one lone electron pair and it needs to be attached to an sp2 hybridized atom. \n- In addition, Hueckel's rule '4n + 2' needs to be true, where 'n' is either zero or any positive integer. \n\nThe corresponding compound structure can be plotted with the ring bonds highlighted in color:","metadata":{}},{"cell_type":"code","source":"mol_num <- 2\nringatoms <- rings(bb_smiles_sdfset[mol_num], type=\"all\")\n\nfig(10, 10)\natomindex <- as.numeric(gsub(\".*_\", \"\", unique(unlist(ringatoms))))\nplot(bb_smiles_sdfset[mol_num], print = FALSE, colbonds = atomindex) ","metadata":{"execution":{"iopub.status.busy":"2024-04-09T15:27:04.653917Z","iopub.execute_input":"2024-04-09T15:27:04.655630Z","iopub.status.idle":"2024-04-09T15:27:04.964422Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ringatoms <- rings(bb_smiles_sdfset[mol_num], type = \"count\", arom = \"TRUE\")\nringatoms","metadata":{"execution":{"iopub.status.busy":"2024-04-09T15:27:04.967168Z","iopub.execute_input":"2024-04-09T15:27:04.968657Z","iopub.status.idle":"2024-04-09T15:27:05.002071Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"If you have ChemmineOB available you can use the regenCoords option to have OpenBabel regenerate the coordinates for the compound. This can sometimes produce better looking plots.","metadata":{}},{"cell_type":"markdown","source":"**Examine the generated features**","metadata":{}},{"cell_type":"code","source":"fig(20, 10)\ncols = c('C', 'H', 'N', 'O', 'Cl', 'S', 'F', 'Br', 'I', 'B', 'Si')\nboxplot(mol_features[, ..cols], col = \"#99BC85\",\n main = \"Atom Frequency in the building blocks SMILES\",\n cex.main = 1.5, cex.axis = 1.5)","metadata":{"execution":{"iopub.status.busy":"2024-04-09T15:27:05.012488Z","iopub.execute_input":"2024-04-09T15:27:05.014165Z","iopub.status.idle":"2024-04-09T15:27:05.349798Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(25, 10)\ncols = c('RNH2', 'R2NH', 'R3N', 'ROPO3', 'ROH', 'RCHO',\n 'RCOR', 'RCOOH', 'RCOOR', 'ROR', 'RCCH', 'RCN',\n 'RINGS', 'AROMATIC')\nboxplot(mol_features[, ..cols], col = \"#99BC85\",\n main = \"Frequency of functional groups, rings, and aromatic rings in the building blocks SMILES\",\n cex.main = 1.5, cex.axis = 1.5) ","metadata":{"execution":{"iopub.status.busy":"2024-04-09T15:27:05.352653Z","iopub.execute_input":"2024-04-09T15:27:05.354087Z","iopub.status.idle":"2024-04-09T15:27:05.704715Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Remove groups absent in all smiles","metadata":{}},{"cell_type":"code","source":"cols <- names(colSums(mol_features[, !c(\"MF\", \"SMILES\")]))[\n colSums(mol_features[, !c(\"MF\", \"SMILES\")]) == 0]\nmol_features <- mol_features[, !cols, with = FALSE]","metadata":{"execution":{"iopub.status.busy":"2024-04-09T15:27:05.709094Z","iopub.execute_input":"2024-04-09T15:27:05.710668Z","iopub.status.idle":"2024-04-09T15:27:05.728325Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(10, 10)\nboxplot(mol_features[, MW], col = \"#99BC85\",\n main = \"Molecular weight of the building blocks SMILES\",\n cex.main = 1.5, cex.axis = 1.5) ","metadata":{"execution":{"iopub.status.busy":"2024-04-09T15:27:05.731045Z","iopub.execute_input":"2024-04-09T15:27:05.732504Z","iopub.status.idle":"2024-04-09T15:27:05.876209Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The features are stored in the dataset: [BELKA: supplementary calcs & data for ML](https://www.kaggle.com/datasets/antoninadolgorukova/belka-supplementary-calcs-and-data-for-ml/data).","metadata":{}},{"cell_type":"code","source":"fwrite(mol_features, \"bb_mol_phys_chem_features.csv\")","metadata":{"execution":{"iopub.status.busy":"2024-04-09T15:27:05.879078Z","iopub.execute_input":"2024-04-09T15:27:05.880692Z","iopub.status.idle":"2024-04-09T15:27:05.896644Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Similarity Comparisons and Searching","metadata":{}},{"cell_type":"markdown","source":"The function `sdf2ap` computes atom pair descriptors for one or many compounds (Carhart, Smith, and Venkataraghavan 1985; Chen and Reynolds 2002). It returns a searchable atom pair database stored in a container of class `APset`, which can be used for structural similarity searching and clustering. As similarity measure, the Tanimoto coefficient or related coefficients can be used.","metadata":{}},{"cell_type":"code","source":"apset <- sdf2ap(bb_smiles_sdfset)\nview(apset[1:3]) ","metadata":{"execution":{"iopub.status.busy":"2024-04-09T15:27:05.899372Z","iopub.execute_input":"2024-04-09T15:27:05.900794Z","iopub.status.idle":"2024-04-09T15:27:06.133698Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# - Duplicates","metadata":{}},{"cell_type":"markdown","source":"The `cmp.duplicated` function can be used to quickly identify very similar compounds in atom pair sets, which will be frequently, **but not necessarily**, identical compounds.","metadata":{}},{"cell_type":"code","source":"any(cmp.duplicated(apset, type = 1))","metadata":{"execution":{"iopub.status.busy":"2024-04-09T15:27:06.136586Z","iopub.execute_input":"2024-04-09T15:27:06.138096Z","iopub.status.idle":"2024-04-09T15:27:06.655241Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dup <- cmp.duplicated(apset)\ncat(\"Found\", length(cid(apset[dup])), \"duplicated instances:\")\ncid(apset[dup])","metadata":{"execution":{"iopub.status.busy":"2024-04-09T15:27:06.666630Z","iopub.execute_input":"2024-04-09T15:27:06.668175Z","iopub.status.idle":"2024-04-09T15:27:07.193636Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dup_clust_dt <- as.data.table(cmp.duplicated(apset, type = 2))\ncat(\"Duplicates forming a cluster of 3 or more molecules:\")\ndup_clust_dt[CLSZ_100 >= 3 ][order(CLID_100)]","metadata":{"execution":{"iopub.status.busy":"2024-04-09T15:27:07.196331Z","iopub.execute_input":"2024-04-09T15:27:07.197698Z","iopub.status.idle":"2024-04-09T15:27:07.764124Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Plot the structure of the duplicates:","metadata":{}},{"cell_type":"code","source":"fig(25, 15)\nplot(bb_smiles_sdfset[c(\"CMP68\", \"CMP133\", \"CMP210\")], print = FALSE)","metadata":{"execution":{"iopub.status.busy":"2024-04-09T15:27:07.766826Z","iopub.execute_input":"2024-04-09T15:27:07.768234Z","iopub.status.idle":"2024-04-09T15:27:08.188423Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Corresponding SMILES:","metadata":{}},{"cell_type":"code","source":"smiles_dict[c(133, 68, 210)]","metadata":{"execution":{"iopub.status.busy":"2024-04-09T15:27:08.191380Z","iopub.execute_input":"2024-04-09T15:27:08.192899Z","iopub.status.idle":"2024-04-09T15:27:08.214596Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"`C=CCC(NC(=O)OCC1c2ccccc2-c2ccccc21)C(=O)O`: [Fmoc-alpha-allyl-DL-glycine](https://pubchem.ncbi.nlm.nih.gov/compound/4052587) \n`C=CC[C@@H](NC(=O)OCC1c2ccccc2-c2ccccc21)C(=O)O`: [Fmoc-D-Allylglycine](https://pubchem.ncbi.nlm.nih.gov/compound/7021747) \n`C=CC[C@H](NC(=O)OCC1c2ccccc2-c2ccccc21)C(=O)O`: [Fmoc-L-Allylglycine](https://pubchem.ncbi.nlm.nih.gov/compound/2734457)","metadata":{}},{"cell_type":"markdown","source":"
\n 💡 Note: These compounds are not the same. They are different stereoisomers of Fmoc-allylglycine, which differ in the stereochemistry of the alpha-carbon. The symbols \"@@\" and \"@\" in SMILES represent stereochemistry, with \"@@\" indicating \"cis\" configuration and \"@\" indicating \"trans\" configuration in this context. In the compounds names, the \"DL\" designation implies a racemic mixture of both D and L enantiomers, where the stereochemistry at the alpha-carbon is not defined. D/L indicate D or L configuration of the alpha-carbon.\n
","metadata":{}},{"cell_type":"markdown","source":"# - Maximum common substructure (MCS)","metadata":{}},{"cell_type":"markdown","source":"Maximum common substructure (MCS) algorithms rank among the most sensitive and accurate methods for measuring structural similarities among small molecules. The FMCS (flexible common substructure) algorithm provides an improved MCS search method that allows atom and/or bond mismatches in the substructures shared among two small molecules. The resulting FMCSs are often larger than strict MCSs, resulting in the identification of more common features in their source structures, as well as a higher sensitivity in detecting weak similarities among compounds.\n\nThe `fmcs` function computes the MCS/FMCS shared among two compounds, which can be highlighted in their structure with the plotMCS function.","metadata":{}},{"cell_type":"code","source":"msc_test <- fmcs(bb_smiles_sdfset[2], bb_smiles_sdfset[2110])\nmsc_test\ncat(\"MCS shared among two molecules:\")\n\nfig(25, 10)\nplotMCS(msc_test)","metadata":{"execution":{"iopub.status.busy":"2024-04-09T15:27:08.217195Z","iopub.execute_input":"2024-04-09T15:27:08.218588Z","iopub.status.idle":"2024-04-09T15:27:08.532568Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The result can be accessed as follows:","metadata":{}},{"cell_type":"code","source":"msc_test[[\"stats\"]][\"Tanimoto_Coefficient\"]","metadata":{"execution":{"iopub.status.busy":"2024-04-09T15:27:08.535456Z","iopub.execute_input":"2024-04-09T15:27:08.536815Z","iopub.status.idle":"2024-04-09T15:27:08.551960Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n 💡 Note: The two molecules have 2 Maximum Common Substructures (MCSs), with the largest one containing 11 atoms. Tanimoto and Overlap coefficients are measures of how similar the two molecules are based on their MCS\n
\n","metadata":{}},{"cell_type":"markdown","source":"Extract Tanimoto coefficients for pairwise comparisons of all molecules (in Kaggle this takes about 3.5 hours, even with parallelization and computing of only the upper triangle of the matrix, so I did run this code on my PC and here just load the precomputed matrix).","metadata":{}},{"cell_type":"code","source":"calculate_similarity <- function(sdfset, measure) {\n \n num_cores <- detectCores()\n registerDoParallel(cores = detectCores())\n\n n_mols <- length(sdfset)\n sdfset_list <- lapply(seq(n_mols-1), function(i) sdfset[i:n_mols])\n \n res_lst <- foreach(batch = sdfset_list, .packages = c(\"fmcsR\")) %dopar% {\n res <- sapply(2:length(batch), function(i) fmcs(batch[1], batch[i])[[\"stats\"]][measure])\n res <- c(rep(0, n_mols - length(res)), res)\n }\n res <- c(unlist(res_lst), rep(0, n_mols))\n similarity_matrix <- matrix(res, nrow = n_mols, ncol = n_mols, byrow = TRUE) \n similarity_matrix[lower.tri(similarity_matrix)] <- \n t(similarity_matrix)[lower.tri(similarity_matrix)]\n diag(similarity_matrix) <- 1\n return(similarity_matrix)\n}","metadata":{"execution":{"iopub.status.busy":"2024-04-09T15:27:08.554483Z","iopub.execute_input":"2024-04-09T15:27:08.555770Z","iopub.status.idle":"2024-04-09T15:27:08.567841Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# library(doParallel)\n\n# tictoc::tic()\n# tam_sim_fmcs <- calculate_similarity(bb_smiles_sdfset, measure = \"Tanimoto_Coefficient\")\n# rownames(tam_sim_fmcs) <- cid(bb_smiles_sdfset)\n# colnames(tam_sim_fmcs) <- cid(bb_smiles_sdfset)\n# tictoc::toc()\n# ~ 3.5h in kaggle :(\ntam_sim_fmcs <- qread(\"/kaggle/input/belka-supplementary-calcs-and-data-for-ml/features/tam_sim_fmcs_bb_smiles.qs\")\n\nhc <- hclust(as.dist(1-tam_sim_fmcs), method = \"single\")\nhead(tam_sim_fmcs)","metadata":{"execution":{"iopub.status.busy":"2024-04-09T15:27:08.570419Z","iopub.execute_input":"2024-04-09T15:27:08.571726Z","iopub.status.idle":"2024-04-09T15:27:09.919054Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"palette_len = 10\nbreaks <- c( \n seq(min(tam_sim_fmcs), 0, length.out = ceiling(palette_len/2) + 1),\n seq(max(tam_sim_fmcs)/palette_len, max(tam_sim_fmcs), length.out = floor(palette_len/2))\n)\ncolors <- colorRampPalette(c(\"darkblue\", \"white\",\"darkred\"))(palette_len)\n\nfig(30, 25)\nheat <- pheatmap(tam_sim_fmcs, fontsize = 16,\n color = colors,\n breaks = breaks,\n angle_col = 90,\n show_rownames = FALSE,\n show_colnames = FALSE,\n cluster_cols = hc,\n cluster_rows = FALSE,\n main = paste(\"Similarity using Tanimoto coefficient based on FMCS for building blocks molecules\"))","metadata":{"execution":{"iopub.status.busy":"2024-04-09T15:27:09.921578Z","iopub.execute_input":"2024-04-09T15:27:09.922872Z","iopub.status.idle":"2024-04-09T15:27:27.887080Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cat(\"The number of molecules per cluster (at the level shown by the red line)\")\n\ncut_heigh = 0.34\n# table(sort(cutree(heat$tree_col, h = cut_heigh)))\n\nfig(30, 20)\nplot(heat$tree_col)\nabline(h = cut_heigh, col = \"red\", lty = 3, lwd = 1)","metadata":{"execution":{"iopub.status.busy":"2024-04-09T15:27:27.889824Z","iopub.execute_input":"2024-04-09T15:27:27.891336Z","iopub.status.idle":"2024-04-09T15:27:28.901098Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"clust <- cutree(heat$tree_col, h = cut_heigh)\nmol_features <- data.table(\"smiles_id\" = cid(bb_smiles_sdfset),\n \"smiles\" = smiles_dict$smiles,\n \"Tanimoto_FMCS_cl\" = clust)\ntable(mol_features$Tanimoto_FMCS_cl)[table(mol_features$Tanimoto_FMCS_cl) > 5]","metadata":{"execution":{"iopub.status.busy":"2024-04-09T15:28:15.680510Z","iopub.execute_input":"2024-04-09T15:28:15.682161Z","iopub.status.idle":"2024-04-09T15:28:15.709692Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"mol_features <- unique(bb_smiles[, !\"set\"])[mol_features, on = \"smiles\"]\nmol_features <- unique(mol_features[, BBs := paste(BB, collapse = \", \"), by = \"smiles\"][, !\"BB\"])","metadata":{"execution":{"iopub.status.busy":"2024-04-09T15:28:18.314623Z","iopub.execute_input":"2024-04-09T15:28:18.316946Z","iopub.status.idle":"2024-04-09T15:28:18.365683Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"mol_features[\n , .(n_clusters = uniqueN(Tanimoto_FMCS_cl)),\n by = c(\"BBs\")]","metadata":{"execution":{"iopub.status.busy":"2024-04-09T15:28:29.224893Z","iopub.execute_input":"2024-04-09T15:28:29.226479Z","iopub.status.idle":"2024-04-09T15:28:29.253108Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can visualize clustering result using dimensional reduction technques.","metadata":{}},{"cell_type":"code","source":"pca <- prcomp(tam_sim_fmcs, center = TRUE, scale = TRUE) \nplt_dt <- data.table(id = rownames(pca$x),\n pca$x[, 1:2])\nplt_dt <- plt_dt[smiles_dict, on = \"id\"][bb_smiles_agg, on = \"smiles\"]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ggplot(plt_dt, aes(x = PC1, y = PC2, color = BB)) +\n geom_point(size = 3) +\n theme_bw(base_size = 22) +\n theme(panel.grid = element_blank()) + #,aspect.ratio = 1\n ggtitle(\"PCA of Tanimoto coefficients based on FMCS for building blocks molecules\")","metadata":{"execution":{"iopub.status.busy":"2024-04-09T16:01:48.221309Z","iopub.execute_input":"2024-04-09T16:01:48.227304Z","iopub.status.idle":"2024-04-09T16:01:49.634659Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ggplot(plt_dt, aes(x = PC1, y = PC2, color = set)) +\n geom_point(size = 3) +\n theme_bw(base_size = 22) +\n theme(panel.grid = element_blank()) + #,aspect.ratio = 1\n ggtitle(\"PCA of Tanimoto coefficients based on FMCS for building blocks molecules\")","metadata":{"execution":{"iopub.status.busy":"2024-04-09T16:02:11.947756Z","iopub.execute_input":"2024-04-09T16:02:11.949536Z","iopub.status.idle":"2024-04-09T16:02:13.083852Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# - Fingerprints-based similarity","metadata":{}},{"cell_type":"markdown","source":"Compound similarity searching can also be done with `FPset` Under method one can choose from several predefined similarity measures including Tanimoto (default), Euclidean, Tversky or Dice. Alternatively, one can pass on custom similarity functions.","metadata":{}},{"cell_type":"markdown","source":"Atom pairs can be converted into binary atom pair fingerprints of fixed length. The function `desc2fp` generates fingerprints from descriptor vectors of variable length such as atom pairs stored in APset or list containers. The `FPset` class stores fingerprints of small molecules in a matrix-like representation where every molecule is encoded as a fingerprint of the same type and length. The obtained fingerprints can be used for structure similarity comparisons, searching and clustering.","metadata":{}},{"cell_type":"code","source":"fpset <- desc2fp(apset)\nfpset\nview(fpset[1:2])","metadata":{"execution":{"iopub.status.busy":"2024-04-09T16:08:07.324539Z","iopub.execute_input":"2024-04-09T16:08:07.326210Z","iopub.status.idle":"2024-04-09T16:08:07.930394Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Similarity between the first two molecules:","metadata":{}},{"cell_type":"code","source":"fpSim(fpset[1], fpset[2], method = \"Tanimoto\") # default method","metadata":{"execution":{"iopub.status.busy":"2024-04-09T16:08:10.117814Z","iopub.execute_input":"2024-04-09T16:08:10.121458Z","iopub.status.idle":"2024-04-09T16:08:10.140367Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tam_sim_ap <- sapply(cid(fpset), function(x) fpSim(x = fpset[x], fpset, sorted = FALSE)) \nhc <- hclust(as.dist(1 - tam_sim_ap), method = \"single\")\nhead(tam_sim_ap)","metadata":{"execution":{"iopub.status.busy":"2024-04-09T16:08:12.016911Z","iopub.execute_input":"2024-04-09T16:08:12.030624Z","iopub.status.idle":"2024-04-09T16:09:13.726197Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"palette_len = 10\nbreaks <- c( \n seq(min(tam_sim_ap), 0, length.out = ceiling(palette_len/2) + 1),\n seq(max(tam_sim_ap)/palette_len, max(tam_sim_ap), length.out = floor(palette_len/2))\n)\ncolors <- colorRampPalette(c(\"darkblue\", \"white\",\"darkred\"))(palette_len)\n\nfig(30, 25)\nheat <- pheatmap(tam_sim_ap, fontsize = 16,\n color = colors,\n breaks = breaks,\n angle_col = 90,\n show_rownames = FALSE,\n show_colnames = FALSE,\n cluster_cols = hc,\n cluster_rows = FALSE,\n main = paste(\"Similarity using Tanimoto coefficient based on atom pair fingerprints for building blocks molecules\"))","metadata":{"execution":{"iopub.status.busy":"2024-04-09T16:09:13.728886Z","iopub.execute_input":"2024-04-09T16:09:13.730334Z","iopub.status.idle":"2024-04-09T16:09:31.479088Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cat(\"The number of molecules per cluster (at the level shown by the red line)\")\n\ncut_heigh = 0.315\n# table(sort(cutree(heat$tree_col, h = cut_heigh)))\n\nfig(30, 20)\nplot(heat$tree_col)\nabline(h = cut_heigh, col = \"red\", lty = 3, lwd = 1)","metadata":{"execution":{"iopub.status.busy":"2024-04-09T16:09:31.483149Z","iopub.execute_input":"2024-04-09T16:09:31.485335Z","iopub.status.idle":"2024-04-09T16:09:32.463566Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"clust <- cutree(heat$tree_col, h = cut_heigh)\nmol_features <- data.table(\"smiles_id\" = cid(bb_smiles_sdfset),\n \"smiles\" = smiles_dict$smiles,\n \"Tanimoto_AP_cl\" = clust)\ntable(mol_features$Tanimoto_AP_cl)[table(mol_features$Tanimoto_AP_cl) > 5]","metadata":{"execution":{"iopub.status.busy":"2024-04-09T16:09:32.466272Z","iopub.execute_input":"2024-04-09T16:09:32.467661Z","iopub.status.idle":"2024-04-09T16:09:32.494366Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"mol_features <- unique(bb_smiles[, !\"set\"])[mol_features, on = \"smiles\"]\nmol_features <- unique(mol_features[, BBs := paste(BB, collapse = \", \"), by = \"smiles\"][, !\"BB\"])","metadata":{"execution":{"iopub.status.busy":"2024-04-09T16:09:32.497035Z","iopub.execute_input":"2024-04-09T16:09:32.498631Z","iopub.status.idle":"2024-04-09T16:09:32.542455Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"mol_features[\n , .(n_clusters = uniqueN(Tanimoto_AP_cl)),\n by = c(\"BBs\")]","metadata":{"execution":{"iopub.status.busy":"2024-04-09T16:09:32.545724Z","iopub.execute_input":"2024-04-09T16:09:32.547167Z","iopub.status.idle":"2024-04-09T16:09:32.575577Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's do again PCA and visualize clustering results.","metadata":{}},{"cell_type":"code","source":"pca <- prcomp(tam_sim_ap, center = TRUE, scale = TRUE) \nplt_dt <- data.table(id = rownames(pca$x),\n pca$x[, 1:2])\nplt_dt <- plt_dt[smiles_dict, on = \"id\"][bb_smiles_agg, on = \"smiles\"]","metadata":{"execution":{"iopub.status.busy":"2024-04-09T16:09:53.980496Z","iopub.execute_input":"2024-04-09T16:09:53.982279Z","iopub.status.idle":"2024-04-09T16:09:59.931638Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ggplot(plt_dt, aes(x = PC1, y = PC2, color = BB)) +\n geom_point(size = 3) +\n theme_bw(base_size = 22) +\n theme(panel.grid = element_blank()) + #,aspect.ratio = 1\n ggtitle(\"PCA of Tanimoto coefficient based on atom pair fingerprints for building blocks molecules\")","metadata":{"execution":{"iopub.status.busy":"2024-04-09T16:10:00.087576Z","iopub.execute_input":"2024-04-09T16:10:00.090353Z","iopub.status.idle":"2024-04-09T16:10:01.208769Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ggplot(plt_dt, aes(x = PC1, y = PC2, color = set)) +\n geom_point(size = 3) +\n theme_bw(base_size = 22) +\n theme(panel.grid = element_blank()) + #,aspect.ratio = 1\n ggtitle(\"PCA of Tanimoto coefficients based on FMCS for building blocks molecules\")","metadata":{"execution":{"iopub.status.busy":"2024-04-09T16:10:06.358382Z","iopub.execute_input":"2024-04-09T16:10:06.360044Z","iopub.status.idle":"2024-04-09T16:10:07.456640Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Storing Atom Pair Fingerprints","metadata":{}},{"cell_type":"markdown","source":"Below is the code to coerce fpset to matrix.","metadata":{}},{"cell_type":"code","source":"tam_apfp_dt <- as.data.table(as.matrix(fpset))\nnames(tam_apfp_dt) <- paste(\"f\", 1:ncol(tam_apfp_dt))\ntam_apfp_dt <- cbind(\"SMILES\" = smiles_dict$smiles, tam_apfp_dt)\nhead(tam_apfp_dt)\nshape(tam_apfp_dt)","metadata":{"execution":{"iopub.status.busy":"2024-04-09T16:10:17.132635Z","iopub.execute_input":"2024-04-09T16:10:17.134193Z","iopub.status.idle":"2024-04-09T16:10:17.655413Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fwrite(tam_apfp_dt, \"bb_mol_apfp_features.csv\")","metadata":{"execution":{"iopub.status.busy":"2024-04-09T16:10:19.921807Z","iopub.execute_input":"2024-04-09T16:10:19.923406Z","iopub.status.idle":"2024-04-09T16:10:19.953033Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Binning Clustering","metadata":{}},{"cell_type":"markdown","source":"Compound libraries can be clustered into discrete similarity groups with the binning clustering function cmp.cluster. The function accepts as input an atom pair (`APset`) or a fingerprint (`FPset`) descriptor database as well as a similarity threshold. The binning clustering result is returned in form of a data frame. Single linkage is used for cluster joining. Because an optimum similarity threshold is often not known, the cmp.cluster function can calculate cluster results for multiple cutoffs in one step with almost the same speed as for a single cutoff. One may force the cmp.cluster function to calculate and store the distance matrix by supplying a file name to the save.distances argument. The generated distance matrix can be loaded and passed on to many other clustering methods available in R, such as the hierarchical clustering function hclust.","metadata":{}},{"cell_type":"code","source":"clusters <- cmp.cluster(db = apset, cutoff = c(0.5, 0.6, 0.7), quiet = TRUE)\nhead(clusters)","metadata":{"execution":{"iopub.status.busy":"2024-04-09T16:10:22.756695Z","iopub.execute_input":"2024-04-09T16:10:22.758282Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cluster.sizestat(clusters, cluster.result = 1)","metadata":{"execution":{"iopub.status.busy":"2024-04-09T11:33:34.784699Z","iopub.execute_input":"2024-04-09T11:33:34.786532Z","iopub.status.idle":"2024-04-09T11:33:34.822048Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cluster.sizestat(clusters, cluster.result = 3)","metadata":{"execution":{"iopub.status.busy":"2024-04-09T11:33:34.825255Z","iopub.execute_input":"2024-04-09T11:33:34.827007Z","iopub.status.idle":"2024-04-09T11:33:34.857435Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n 💡 Note: \n The clustering algorithm produced varying numbers of clusters based on different cutoff values: lower cutoff values result in larger numbers of smaller clusters, while higher cutoff values lead to fewer but larger clusters. Optimal cutoff values often lead to a reasonable balance between the number of clusters and their sizes.\n
","metadata":{}},{"cell_type":"markdown","source":"We can also visualize clustering result using multi-dimensional scaling (MDS).","metadata":{}},{"cell_type":"code","source":"cluster.visualize(apset, clusters, size.cutoff = 1, quiet = TRUE, cluster.result = 1)","metadata":{"execution":{"iopub.status.busy":"2024-04-09T11:33:34.860660Z","iopub.execute_input":"2024-04-09T11:33:34.862434Z","iopub.status.idle":"2024-04-09T11:35:40.950679Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n 💡 Note: It seems that all clustering methods show a very similar picture with 2 relatively large clusters and many small ones. The fact that all BB smiles form consistently across different clustering methods suggests that these clusters may represent meaningful patterns in the data.\n