{"metadata":{"kernelspec":{"name":"ir","display_name":"R","language":"R"},"language_info":{"name":"R","codemirror_mode":"r","pygments_lexer":"r","mimetype":"text/x-r-source","file_extension":".r","version":"4.0.5"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Come on Chemicals! Chemical structure based EDA in R\n\nThis is an attempt to drill down into potential feature engineering around chemicals. Heavily inspired by [this feature engineering notebook](https://www.kaggle.com/code/antoninadolgorukova/op2-feature-engineering).\n\nThe main goal of this notebook is to determine whether or not quantitative measures of chemical structure seem to generalize to the gene perturbation space.\n\nWork by Vendekagon Labs. If this kind of first principles based feature investigation is interesting to you, you might also want to check out the [single cell EDA notebook](https://www.kaggle.com/code/vendekagonlabs/op2-single-cell-eda-10x-multiome).","metadata":{}},{"cell_type":"code","source":"install.packages(\"rcdk\", verbose = FALSE)\nlibrary(\"arrow\")\nlibrary(\"data.table\")\nlibrary(\"dendextend\")\nlibrary(\"qs\")\nlibrary(\"rcdk\")\nlibrary(\"fingerprint\") \nlibrary(\"plotly\")\nlibrary(\"umap\")\nlibrary(\"dbscan\")","metadata":{"_uuid":"051d70d956493feee0c6d64651c6a088724dca2a","_execution_state":"idle","execution":{"iopub.status.busy":"2023-09-24T16:27:02.093283Z","iopub.execute_input":"2023-09-24T16:27:02.095431Z","iopub.status.idle":"2023-09-24T16:27:13.966551Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"list.files(path = \"/kaggle/input/open-problems-single-cell-perturbations/\")","metadata":{"execution":{"iopub.status.busy":"2023-09-24T16:27:13.969587Z","iopub.execute_input":"2023-09-24T16:27:13.970925Z","iopub.status.idle":"2023-09-24T16:27:13.995637Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Reading in training data\n\nReading in a parquet using the arrow library and converting to a `data.table`.","metadata":{}},{"cell_type":"code","source":"de_train <- read_parquet('/kaggle/input/open-problems-single-cell-perturbations/de_train.parquet')\nde_train <- setDT(de_train)\nhead(de_train)","metadata":{"execution":{"iopub.status.busy":"2023-09-24T16:27:13.998307Z","iopub.execute_input":"2023-09-24T16:27:13.999738Z","iopub.status.idle":"2023-09-24T16:27:21.142288Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Splitting into gene expression change based features and the small molecule names for later use.","metadata":{}},{"cell_type":"code","source":"feats = de_train[,6:length(colnames(de_train))]\nchems = de_train$sm_name","metadata":{"execution":{"iopub.status.busy":"2023-09-24T16:27:21.145153Z","iopub.execute_input":"2023-09-24T16:27:21.146504Z","iopub.status.idle":"2023-09-24T16:27:21.178102Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Investigating Chemical Structure\n\nThis short section follows the other notebook fairly heavily. The goal here is to look at a variety of ways to encode the molecular structure quantitatively, as well as to characterize similarity.","metadata":{}},{"cell_type":"code","source":"mols <- parse.smiles(unique(de_train$SMILES))\nfps <- lapply(mols, get.fingerprint, type='circular')\nnames(fps) <- unique(de_train$sm_name)\nhead(fps)","metadata":{"execution":{"iopub.status.busy":"2023-09-24T16:27:21.180248Z","iopub.execute_input":"2023-09-24T16:27:21.181370Z","iopub.status.idle":"2023-09-24T16:27:21.450780Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Producing `finger_features`, which is a binary encoding of the chemical structure of each compound.","metadata":{}},{"cell_type":"code","source":"finger_features <- cbind(data.table(\"sm_name\" = names(fps)),\n                        fp.to.matrix(fps))\nhead(finger_features)","metadata":{"execution":{"iopub.status.busy":"2023-09-24T16:27:21.454019Z","iopub.execute_input":"2023-09-24T16:27:21.455886Z","iopub.status.idle":"2023-09-24T16:27:21.867018Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ff_tab <- setDT(finger_features)\nfwrite(ff_tab, \"/kaggle/working/finger_features.csv\")","metadata":{"execution":{"iopub.status.busy":"2023-09-24T16:27:21.869524Z","iopub.execute_input":"2023-09-24T16:27:21.871064Z","iopub.status.idle":"2023-09-24T16:27:21.887107Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Hierarchical Clustering of Chemical Structure\n\nStill following the [other notebook](https://www.kaggle.com/code/antoninadolgorukova/op2-feature-engineering),\nwe cluster molecules based on their Tanimoto similarity metric. Then use a heuristic method to derive clusters of molecules.","metadata":{}},{"cell_type":"code","source":"options(repr.plot.width = 25, repr.plot.height = 10)\nfp.sim <- fingerprint::fp.sim.matrix(fps, method='tanimoto')\nfp.dist <- 1 - fp.sim\ncls <- hclust(as.dist(fp.dist))\ncls$labels <- names(fps)\nplot(cls, main = 'Clustering of the compounds based on their molecular fingerprints and the Tanimoto similarity metric')","metadata":{"execution":{"iopub.status.busy":"2023-09-24T16:27:21.889751Z","iopub.execute_input":"2023-09-24T16:27:21.891028Z","iopub.status.idle":"2023-09-24T16:27:22.208406Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#------------\n# this code block verbatim from:\n# https://www.kaggle.com/code/antoninadolgorukova/op2-feature-engineering\n# -----------\nheights <- cls$height\nn_clusters <- length(heights) + 1  # Number of clusters including individual data points\nwss <- numeric(n_clusters)\n\n# Calculate within-cluster sum of squares (WSS) for different levels\nfor (k in 1:n_clusters) {\n  cut_tree <- cutree(cls, k = k)\n  cluster_heights <- heights[cut_tree]  # Heights of clusters at level k\n  wss[k] <- sum((cluster_heights - mean(cluster_heights))^2)  # WSS calculation\n}\n\n# Find the \"Elbow Point\" using a heuristic\nelbow_point <- 1\nfor (k in 2:(n_clusters - 1)) {\n  if (wss[k] < wss[k - 1] && wss[k] < wss[k + 1]) {\n    elbow_point <- k\n    break\n  }\n}\n\n# Visualize the \"Elbow Point\" on the dendrogram\ndend <- as.dendrogram(cls)\ndend %>% \n    set(\"labels_col\", rainbow(elbow_point), k = elbow_point) %>% # change color\n    set(\"branches_k_color\", k = elbow_point) %>% \n  plot(main = 'Clustering of the compounds based on Elbow point (k = 11)')\nabline(h = heights[elbow_point], col = \"red\", lty = 2)","metadata":{"execution":{"iopub.status.busy":"2023-09-24T16:31:22.107055Z","iopub.execute_input":"2023-09-24T16:31:22.113573Z","iopub.status.idle":"2023-09-24T16:31:22.695589Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"clusters <- cutree(cls, k = elbow_point)\nfinger_features[, cluster:= clusters]\n\nfor (cl in unique(finger_features$cluster)) {\n    \n    cat(\"\\ncluster №\", cl, \"\\n\")\n    print(finger_features[cluster == cl, sm_name])\n}","metadata":{"execution":{"iopub.status.busy":"2023-09-24T16:31:25.453299Z","iopub.execute_input":"2023-09-24T16:31:25.456728Z","iopub.status.idle":"2023-09-24T16:31:25.517063Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Additional Cluster Investigation\n\nWe explore the similarity space as well as the binary encoding space for the chemical structures with UMAP projections, and examine the placement of the clusters in the space. We also explore HDBSCAN as an alternative hiearchical clustering algorithm, applied in the UMAP space. We get some fairly similar clustering results.","metadata":{}},{"cell_type":"code","source":"library('umap')\nfp.umap = umap(fp.sim, n_components = 3)\nlayout <- fp.umap[[\"layout\"]]\nlayout <- data.frame(layout)\nclusters <- data.frame(clusters)\nlayout <- cbind(layout, clusters$clusters)\nlayout <- cbind(layout, rownames(clusters))\ncolnames(layout)[colnames(layout) == \"clusters$clusters\"] =\"cluster\"\ncolnames(layout)[colnames(layout) == \"rownames(clusters)\"] =\"molecule\"\nlayout$cluster_name <- as.factor(layout$cluster)\nhead(layout)","metadata":{"execution":{"iopub.status.busy":"2023-09-24T16:27:22.781786Z","iopub.execute_input":"2023-09-24T16:27:22.783088Z","iopub.status.idle":"2023-09-24T16:27:23.335272Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We have a lot of clusters! So as a quick hack, we stitch together two color palettes to accommodate our number of clusters.","metadata":{}},{"cell_type":"code","source":"library(\"RColorBrewer\")\ndisplay.brewer.all(type = 'qual')\npal1 <- brewer.pal(12, 'Set3')\npal2 <- brewer.pal(12, 'Paired')\nbigpal <- c(pal1, pal2)\nbigpal","metadata":{"execution":{"iopub.status.busy":"2023-09-24T16:29:25.983675Z","iopub.execute_input":"2023-09-24T16:29:25.985101Z","iopub.status.idle":"2023-09-24T16:29:26.220168Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_ly(layout, x = ~X1, y = ~X2, z = ~X3, type=\"scatter3d\", mode='markers',\n        color= ~cluster_name, text = ~cluster_name, colors = bigpal)","metadata":{"execution":{"iopub.status.busy":"2023-09-24T16:28:05.122852Z","iopub.execute_input":"2023-09-24T16:28:05.126293Z","iopub.status.idle":"2023-09-24T16:28:05.934160Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ff.umap = umap(finger_features[,2:length(colnames(finger_features))] , n_components = 3)\nlayout <- ff.umap[[\"layout\"]]\nlayout <- data.frame(layout)\nclusters <- data.frame(clusters)\nlayout <- cbind(layout, clusters$clusters)\nlayout <- cbind(layout, rownames(clusters))\ncolnames(layout)[colnames(layout) == \"clusters$clusters\"] =\"cluster\"\ncolnames(layout)[colnames(layout) == \"rownames(clusters)\"] =\"molecule\"\nlayout$cluster_name <- as.factor(layout$cluster)\nhead(layout)","metadata":{"execution":{"iopub.status.busy":"2023-09-24T16:27:24.605317Z","iopub.execute_input":"2023-09-24T16:27:24.606355Z","iopub.status.idle":"2023-09-24T16:27:25.228815Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig3 <- plot_ly(layout, x = ~X1, y = ~X2, z = ~X3, type=\"scatter3d\", mode='markers',\n                colors = bigpal, color= ~cluster_name)\nfig3","metadata":{"execution":{"iopub.status.busy":"2023-09-24T16:27:25.231237Z","iopub.execute_input":"2023-09-24T16:27:25.232503Z","iopub.status.idle":"2023-09-24T16:27:26.192889Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We then visualize the chemical names directly in plotly 3d. _Note_: for this plot to be useful, you\nshould probably zoom in.","metadata":{}},{"cell_type":"code","source":"cl <- hdbscan(layout[,1:3], minPts=3)\nlayout$hbscan_cluster <- as.factor(cl$cluster)\nplot_ly(layout, x = ~X1, y = ~X2, z = ~X3, type=\"scatter3d\", mode='text',\n        color= ~hbscan_cluster, text = ~molecule, colors = bigpal)","metadata":{"execution":{"iopub.status.busy":"2023-09-24T16:27:26.194770Z","iopub.execute_input":"2023-09-24T16:27:26.195767Z","iopub.status.idle":"2023-09-24T16:27:27.827829Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"just_cls <- layout[,c('molecule', 'hbscan_cluster', 'cluster_name')]\nhead(just_cls)","metadata":{"execution":{"iopub.status.busy":"2023-09-24T16:27:27.830306Z","iopub.execute_input":"2023-09-24T16:27:27.831709Z","iopub.status.idle":"2023-09-24T16:27:27.852203Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sub <- just_cls[just_cls$hbscan_cluster == 18,]\nsub","metadata":{"execution":{"iopub.status.busy":"2023-09-24T16:27:27.854434Z","iopub.execute_input":"2023-09-24T16:27:27.855506Z","iopub.status.idle":"2023-09-24T16:27:27.878587Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Qualitative breakdown:\n\n- Clotrimazole <- anti-fungal\n- Crizotinib <- kinase inhibitor\n- Palbociclib <- CDK 4/6 (kinase) inhibitor\n- LDN 193189 <- ALK2 and ALK3 (kinase) inhibitor\n- CGP 60474 <- CDK (kinase) inhibitor\n- RN-486 <- Btk (kinase) inhibitor\n- Dovitinib <- PTK (kinase) inhibitor\n- K-02288 <- \"K 02288 is a potent and selective inhibitor of type I bone morphogenic protein (BMP) receptors\"\n- BMS-536924 <-  \"ATP-competitive IGF-1R/IR inhibitor\"\n- Dactolisib <- dual ATP-competitive PI3K (kinase) and mTOR inhibitor\n- CHIR-99021 <- GSK-3α and GSK-3β (kinase) inhibitor \n- Tipifarnib <- farnesyltransferase inhibitor (blocks Ras binding to membrane)\n\nAlso, `-nib` suffix denotes small molecule inhibitor.\n\nLooking at the outlier:\n\n> Clotrimazole exerts its action primarily by damaging the permeability barrier in the fungal cytoplasmic membrane.[6] Clotrimazole thereby inhibits the biosynthesis of ergosterol in a concentration-dependent manner by inhibiting the demethylation of 14 alpha lanosterol.\n\n🤷","metadata":{}},{"cell_type":"code","source":"sub <- just_cls[just_cls$hbscan_cluster == 6,]\nsub","metadata":{"execution":{"iopub.status.busy":"2023-09-24T16:27:27.881413Z","iopub.execute_input":"2023-09-24T16:27:27.882663Z","iopub.status.idle":"2023-09-24T16:27:27.901760Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- Lamivudine <- reverse transcriptase inhibitor (so antiviral, re: retrovirus)\n- Decitabine <- hypomethylation agent (specifically, kills abnormal cells in blood marrow)\n- Mubritinib <- HER2 inhibitor (cancer therapeutic)\n- Azacitidine <- chemo (anti-metabolite)","metadata":{}},{"cell_type":"code","source":"sub <- just_cls[just_cls$hbscan_cluster == 7,]\nsub","metadata":{"execution":{"iopub.status.busy":"2023-09-24T16:27:27.905900Z","iopub.execute_input":"2023-09-24T16:27:27.907682Z","iopub.status.idle":"2023-09-24T16:27:27.930032Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"- Dabrafenib <- Kinase inhibitor\n- Oxybenzone <- suncreen component, controversial evidence re: environment harm, allergy, etc.\n- TGX 221 <- Kinase inhibitor\n- SB525334 <- TGFβ receptor I (ALK5) inhibitor","metadata":{}},{"cell_type":"code","source":"feat.umap = umap(feats , n_components = 3)\nfproj <- feat.umap[[\"layout\"]]\nfproj <- data.frame(fproj)\nhead(fproj)","metadata":{"execution":{"iopub.status.busy":"2023-09-24T16:27:27.932462Z","iopub.execute_input":"2023-09-24T16:27:27.933624Z","iopub.status.idle":"2023-09-24T16:27:39.100693Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Back to Gene Expression Space\n\nSo there are maybe some interesting patterns in the chemical space, but it doesn't seem like we can consistently relate them back to drugs and their mechanism of action. Some clusters seem to have a theme, others don't. What happens when we go back into gene expression space with our cluster information?","metadata":{}},{"cell_type":"code","source":"fproj$chems = chems\nplot_ly(fproj, x = ~X1, y = ~X2, z = ~X3, type=\"scatter3d\",\n        mode = 'markers', text = ~chems)","metadata":{"execution":{"iopub.status.busy":"2023-09-24T16:27:39.103087Z","iopub.execute_input":"2023-09-24T16:27:39.104220Z","iopub.status.idle":"2023-09-24T16:27:39.666511Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fproj$molecule <- fproj$chems\nmerged <- merge(fproj, just_cls, by='molecule')\nhead(merged)","metadata":{"execution":{"iopub.status.busy":"2023-09-24T16:27:39.669864Z","iopub.execute_input":"2023-09-24T16:27:39.671742Z","iopub.status.idle":"2023-09-24T16:27:39.706018Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_ly(merged, x = ~X1, y = ~X2, z = ~X3, type = \"scatter3d\", colors = bigpal,\n        mode='markers', text = ~molecule, color = ~hbscan_cluster)","metadata":{"execution":{"iopub.status.busy":"2023-09-24T16:27:39.956658Z","iopub.execute_input":"2023-09-24T16:27:39.958297Z","iopub.status.idle":"2023-09-24T16:27:40.729758Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Exploring in the space above, I see a small cluster of chemicals which includes Scriptaid and Vorinostat. This is interesting, because Scriptaid led to the discovery of the structurally similar Vorinostat according to the Scriptaid [wikipedia article](https://en.wikipedia.org/wiki/Scriptaid). And these are very, very close in feature space with gene expression perturbation. But they're in different clusters in our chemical encoding space with either HDBSCAN or standard hierarchical clustering.","metadata":{}},{"cell_type":"code","source":"new_sub <- merged[merged$molecule %in% c(\"Scriptaid\", \"Vorinostat\"),]\nnew_sub","metadata":{"execution":{"iopub.status.busy":"2023-09-24T16:27:40.732038Z","iopub.execute_input":"2023-09-24T16:27:40.733160Z","iopub.status.idle":"2023-09-24T16:27:40.758636Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_ly(new_sub, x = ~X1, y = ~X2, z = ~X3, type=\"scatter3d\", colors = bigpal,\n        mode = 'markers', text = ~molecule, color = 'cluster_name')","metadata":{"execution":{"iopub.status.busy":"2023-09-24T16:27:40.760812Z","iopub.execute_input":"2023-09-24T16:27:40.761907Z","iopub.status.idle":"2023-09-24T16:27:41.402865Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Average Gene Expression Perturbation Embeddings\n\nFor one last check, we average the gene expression perturbation UMAP embedding vectors.","metadata":{}},{"cell_type":"code","source":"to_merge <- merged[,c('cluster_name', 'molecule')]\nnames(to_merge) <- c('cluster_name', 'molecule')\nhead(to_merge)","metadata":{"execution":{"iopub.status.busy":"2023-09-24T16:27:41.405369Z","iopub.execute_input":"2023-09-24T16:27:41.406615Z","iopub.status.idle":"2023-09-24T16:27:41.427374Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"avg_proj <- aggregate(merged[,c('X1', 'X2', 'X3')], by = list(merged$molecule), FUN=mean)\nnames(avg_proj) <- c('molecule', 'X1', 'X2', 'X3')\navg_proj <- merge(avg_proj, unique(to_merge), by='molecule')\nplot_ly(avg_proj, x = ~X1, y = ~X2, z = ~X3, type = \"scatter3d\", mode='markers',\n        colors = bigpal, text = ~molecule, color= ~cluster_name)","metadata":{"execution":{"iopub.status.busy":"2023-09-24T16:27:41.430272Z","iopub.execute_input":"2023-09-24T16:27:41.431771Z","iopub.status.idle":"2023-09-24T16:27:42.243204Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Conclusions\n\nVisualizing as per above, and with some other spot checking, it doesn't really look like the clusters provide us much structure in the gene perturbation space. It probably makes sense to look elsewhere for feature engineering or chemical grouping, perhaps in embeddings from a deep learning model?","metadata":{}}]}