Calculating activity level for all train building blocks SMILES\n
Exploring do the position of a building block within a molecule plays a role in determining its binding activity\n
Trying to find what characterizes a highly active building blocks\n
\n In this notebook I reproduce approach used in this article with some additional considerations. \n
","metadata":{}},{"cell_type":"code","source":"suppressPackageStartupMessages({\n library(arrow)\n library(duckdb)\n library(data.table)\n library(qs)\n library(ggplot2)\n library(foreach)\n library(doParallel)\n})","metadata":{"_uuid":"051d70d956493feee0c6d64651c6a088724dca2a","_execution_state":"idle","execution":{"iopub.status.busy":"2024-06-29T04:52:25.025676Z","iopub.execute_input":"2024-06-29T04:52:25.027817Z","iopub.status.idle":"2024-06-29T04:52:26.696127Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig <- function(width, heigth) {\n options(repr.plot.width = width, repr.plot.height = heigth)\n}","metadata":{"execution":{"iopub.status.busy":"2024-06-29T04:52:26.698840Z","iopub.execute_input":"2024-06-29T04:52:26.748158Z","iopub.status.idle":"2024-06-29T04:52:26.762550Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Load train data and connect to the DuckDB**","metadata":{}},{"cell_type":"code","source":"train_meta <- open_dataset('/kaggle/input/leash-BELKA/train.parquet')\n\ncon <- dbConnect(duckdb::duckdb())\ntrain <- arrow::to_duckdb(train_meta, table_name = \"train\", con = con)\ntrain","metadata":{"execution":{"iopub.status.busy":"2024-06-29T04:52:26.765311Z","iopub.execute_input":"2024-06-29T04:52:26.766830Z","iopub.status.idle":"2024-06-29T04:52:28.478863Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"proteins <- c(\"BRD4\", \"HSA\", \"sEH\")","metadata":{"execution":{"iopub.status.busy":"2024-06-29T04:52:28.483306Z","iopub.execute_input":"2024-06-29T04:52:28.485342Z","iopub.status.idle":"2024-06-29T04:52:28.503525Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Extract smiles for all building blocks**","metadata":{}},{"cell_type":"code","source":"bb_smiles <- unique(rbindlist(lapply(1:3, function(i) {\n as.data.table(\n dbGetQuery(con, paste0(\n \"SELECT DISTINCT(buildingblock\", i, \"_smiles) AS smiles FROM train\"\n ))\n )[, BB := paste0(\"buildingblock\", i, \"_smiles\")]\n})))\n\nbbs_cols <- bb_smiles[, unique(BB)]","metadata":{"execution":{"iopub.status.busy":"2024-06-29T04:52:28.506912Z","iopub.execute_input":"2024-06-29T04:52:28.508574Z","iopub.status.idle":"2024-06-29T04:53:16.343093Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Building blocks (BB) activity level","metadata":{}},{"cell_type":"markdown","source":"To calculate P(active) scores, for each BB, I count the number of molecules that bind to the protein when it occurs at a given position (num_a) and the total number of molecules with that BB at that position (total_mols). Next, I calculate P(active) - the measure of activity - as num_a/total_mols, convert to percentage (x100) and round to make more readable. This took about an hour on my PC, so to save time here I load the precalculated file.","metadata":{}},{"cell_type":"code","source":"calc_p_active <- function(con, bbs_col, smiles) {\n \n num_a <- dbGetQuery(con, paste0(\"\n SELECT protein_name, SUM(binds) AS num_active\n FROM train\n WHERE \", bbs_col, \" = '\", smiles, \"'\n GROUP BY protein_name\n \"\n ))\n \n total_mols <- dbGetQuery(con, paste0(\"\n SELECT COUNT(*) AS n_mols\n FROM train\n WHERE \", bbs_col, \" = '\", smiles, \"'\n \"\n ))\n \n res_dt <- as.data.table(num_a)\n res_dt <- res_dt[, .(protein_name, BB = bbs_col,\n SMILES = smiles, num_active,\n n_mols = total_mols$n_mols/3)]\n \n return(res_dt)\n}\n# tictoc::tic()\n# bb_summary <- rbindlist(\n# lapply(bbs_cols, function(bbs_col) {\n \n# res <- rbindlist(\n# lapply(bb_smiles[BB == bbs_col]$smiles[1:2], function(smiles) {\n# calc_p_active(con, bbs_col, smiles)\n# })\n# )\n# })\n# )\n# tictoc::toc()\n# 4297.917 sec elapsed","metadata":{"execution":{"iopub.status.busy":"2024-06-29T04:53:16.348719Z","iopub.execute_input":"2024-06-29T04:53:16.351412Z","iopub.status.idle":"2024-06-29T04:53:16.367072Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# qsave(bb_summary, \"gen_data/features/bb_summary.qs\")\nbb_summary <- qread(\"/kaggle/input/belka-supplementary-calcs-and-data-for-ml/smiles/bb_summary.qs\")\n\nbb_summary[, p_active := round(num_active/n_mols * 100, 2)]\nhead(bb_summary)","metadata":{"execution":{"iopub.status.busy":"2024-06-29T04:53:16.369803Z","iopub.execute_input":"2024-06-29T04:53:16.371309Z","iopub.status.idle":"2024-06-29T04:53:16.429759Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# The role of the building block (BB) position","metadata":{}},{"cell_type":"markdown","source":"Reshape the bb_summary data table to wide format to see percentage of active molecules (P(active)) respective to the three proteins:","metadata":{}},{"cell_type":"code","source":"bb_wide <- dcast(bb_summary, BB + SMILES ~ protein_name, value.var = \"p_active\")\nhead(bb_wide)","metadata":{"execution":{"iopub.status.busy":"2024-06-29T04:53:16.432658Z","iopub.execute_input":"2024-06-29T04:53:16.434380Z","iopub.status.idle":"2024-06-29T04:53:16.478265Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Range of the P(active) score:**","metadata":{}},{"cell_type":"code","source":"bb_summary[, range(p_active)]","metadata":{"execution":{"iopub.status.busy":"2024-06-29T04:53:16.480969Z","iopub.execute_input":"2024-06-29T04:53:16.482584Z","iopub.status.idle":"2024-06-29T04:53:16.501565Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Are there BBs, with which no single molecule binds to all proteins?**","metadata":{}},{"cell_type":"code","source":"bb_wide[rowSums(bb_wide[, ..proteins]) == 0]","metadata":{"execution":{"iopub.status.busy":"2024-06-29T04:53:16.504196Z","iopub.execute_input":"2024-06-29T04:53:16.505617Z","iopub.status.idle":"2024-06-29T04:53:16.532397Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"bb_summary[SMILES == \"NCc1ccccn1\"][order(BB)]","metadata":{"execution":{"iopub.status.busy":"2024-06-29T04:53:16.535130Z","iopub.execute_input":"2024-06-29T04:53:16.536744Z","iopub.status.idle":"2024-06-29T04:53:16.579754Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"bb_summary[SMILES == \"O=C(Nc1c(Cl)cc(Cl)nc1C(=O)O)OCC1c2ccccc2-c2ccccc21\"][order(BB)]","metadata":{"execution":{"iopub.status.busy":"2024-06-29T04:53:16.582378Z","iopub.execute_input":"2024-06-29T04:53:16.583829Z","iopub.status.idle":"2024-06-29T04:53:16.612975Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"bb_summary[SMILES == \"O=C(Nc1cc(Cl)ncc1C(=O)O)OCC1c2ccccc2-c2ccccc21\"][order(BB)]","metadata":{"execution":{"iopub.status.busy":"2024-06-29T04:53:16.615645Z","iopub.execute_input":"2024-06-29T04:53:16.617044Z","iopub.status.idle":"2024-06-29T04:53:16.644014Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n 💡 Note: So with 'NCc1ccccn1' in BB №3 none of the 266 molecules binds to any protein. Also, when this BB is in 2nd position (BB №2, there are 235,364 such molecules), some molecules binds to proteins, but the proportion is very small (up to 0.5%). For the rest two BBs we see p_active = 0, because of rounding, but anyway, only negligible amount of molecules with them binds to any of the proteins (2-15 of 362k).\n
","metadata":{}},{"cell_type":"markdown","source":"**Are there BBs with position-specific activity level?**","metadata":{}},{"cell_type":"markdown","source":"To answer this, I subset BBs present in both 2nd and third positions (buildingblock2_smiles and buildingblock3_smiles columns). Then I sort the result and calculate percent difference in activity score when each BB (SMILES) localed in 2nd vs 3rd position (BB).","metadata":{}},{"cell_type":"code","source":"common_bbs <- intersect(bb_wide[BB == \"buildingblock2_smiles\", SMILES],\n bb_wide[BB == \"buildingblock3_smiles\", SMILES])\nlength(common_bbs)\nact_diffs <- bb_wide[SMILES %in% common_bbs][order(SMILES, BB)]\nact_diffs <- act_diffs[, lapply(.SD, diff), by = c(\"SMILES\"), .SD = proteins]\nhead(act_diffs)","metadata":{"execution":{"iopub.status.busy":"2024-06-29T04:53:16.646634Z","iopub.execute_input":"2024-06-29T04:53:16.648057Z","iopub.status.idle":"2024-06-29T04:53:16.718242Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Negative values indicate that the percentage of binding molecules (P(active)) is smaller when the BB is located in the 3rd position compared to the 2nd position. Let's check what if there is a more successful position for all BB and proteins and the ranges of the differences:","metadata":{}},{"cell_type":"code","source":"lapply(proteins, function(p) {\n data.table(protein = p,\n n_3rd_higher = act_diffs[get(p) > 0, .N],\n n_3rd_smaller = act_diffs[get(p) < 0, .N],\n n_3rd_no_diff = act_diffs[get(p) == 0, .N]\n )\n})","metadata":{"execution":{"iopub.status.busy":"2024-06-29T04:53:16.721003Z","iopub.execute_input":"2024-06-29T04:53:16.722561Z","iopub.status.idle":"2024-06-29T04:53:16.779089Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"act_diffs[, lapply(.SD, range), .SDcols = proteins]","metadata":{"execution":{"iopub.status.busy":"2024-06-29T04:53:16.781692Z","iopub.execute_input":"2024-06-29T04:53:16.783125Z","iopub.status.idle":"2024-06-29T04:53:16.806209Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n 💡 Note: \n
\n
For each protein (BRD4, HSA, sEH), I've calculated the counts of instances where the activity score is higher, smaller, or not different when the BB is in the 3rd position compared to the 2nd position, e.g. for BRD4, there are 255 instances where the activity score is higher when a BB is in the 3rd position, 405 instances where it is smaller, and 31 instances where there is no difference.\n
So we see that the binding status can vary depending on both the target (protein) and the position of the BB.\n
Obviolsy, the surrounding context within the molecule influenced the activity scores I calculated for each BB, but the train data lacks molecules with the same building block (BB) in the first position and two identical BBs but in different order in the 2nd and 3rd positions (like ABC and ACB) to explore this.\n
\n
","metadata":{}},{"cell_type":"markdown","source":"**How the activity of molecules with the same BB x position combination differ between proteins?**","metadata":{}},{"cell_type":"code","source":"bb_wide[, max_diff := pmax(\n abs(BRD4 - HSA), \n abs(BRD4 - sEH), \n abs(HSA - sEH)\n)]\nbb_wide[order(-max_diff)][1:10]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"
\n 💡 Note: \n
\n
I filtered 10 rows with the maximal difference in any combination of the three protein columns\n
While there are molecules with specific BB x positon combination which do not bind one or two of the proteins, other molecules with the same BB x positon combination bind other proteins (indicates that the relationship between BB x position combinations and protein binding is protein-specific.)\n
\n
","metadata":{}},{"cell_type":"markdown","source":"**Minimal and maximal P(active) score for each protein:**","metadata":{}},{"cell_type":"code","source":"bb_summary[, .(min_p_active = min(p_active),\n max_p_active = max(p_active)),\n by = \"protein_name\"]","metadata":{"execution":{"iopub.status.busy":"2024-06-29T05:00:06.343536Z","iopub.execute_input":"2024-06-29T05:00:06.345268Z","iopub.status.idle":"2024-06-29T05:00:06.377811Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Maximal P(active) score for each protein by BB position:**","metadata":{}},{"cell_type":"code","source":"max_active <- bb_summary[, .(max_p_active = max(p_active)),\n by = c(\"protein_name\", \"BB\")]\n\ndcast(max_active, ... ~ protein_name, value.var = \"max_p_active\", fill = 0)","metadata":{"execution":{"iopub.status.busy":"2024-06-29T05:05:41.346749Z","iopub.execute_input":"2024-06-29T05:05:41.348498Z","iopub.status.idle":"2024-06-29T05:05:41.381385Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Interesting is 71.95% for BB1 and sEH, since BB1 has the smallest effect on overall compound activity as the position closest to the DNA tag (obstruct interactions with the target)","metadata":{}},{"cell_type":"code","source":"bb_summary[p_active > 70]","metadata":{"execution":{"iopub.status.busy":"2024-06-29T05:06:23.858376Z","iopub.execute_input":"2024-06-29T05:06:23.860996Z","iopub.status.idle":"2024-06-29T05:06:23.896372Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(25, 8)\n\nggplot(bb_summary, aes(x = p_active)) + \n geom_histogram(bins = 30, color = \"black\", fill = \"steelblue\") +\n facet_wrap(~ protein_name, scales = \"free\") +\n ggtitle(\"Distribution of P(active) across different proteins\") +\n theme_bw(base_size = 22) +\n theme(legend.position = \"bottom\",\n panel.grid.major = element_line(linewidth = 0.2),\n panel.grid.minor = element_line(linewidth = 0.2),\n panel.border = element_rect(linewidth = 0.5),\n strip.background = element_blank()\n )","metadata":{"execution":{"iopub.status.busy":"2024-06-29T05:08:08.731799Z","iopub.execute_input":"2024-06-29T05:08:08.733613Z","iopub.status.idle":"2024-06-29T05:08:09.883625Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Set P(active) intervals","metadata":{}},{"cell_type":"code","source":"breaks <- seq(0, 100, by = 10)\nbb_summary[, P_act_int := cut(\n p_active, breaks = breaks, dig.lab = 2, include.lowest = TRUE\n )\n]","metadata":{"execution":{"iopub.status.busy":"2024-06-29T05:09:40.242896Z","iopub.execute_input":"2024-06-29T05:09:40.245111Z","iopub.status.idle":"2024-06-29T05:09:40.266016Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(15, 8)\n\nggplot(bb_summary, aes(x = P_act_int)) +\n geom_bar(stat = \"count\", color = \"black\", fill = \"steelblue\") +\n geom_text(stat = \"count\", aes(label = after_stat(count)), vjust = -0.5, size = 6) +\n ggtitle(\"Numbers of BB x position combinations at different P(active) intervals\") +\n xlab(\"P(active) interval\") +\n theme_bw(base_size = 22) +\n theme(legend.position = \"bottom\",\n panel.grid.major = element_line(linewidth = 0.2),\n panel.grid.minor = element_line(linewidth = 0.2),\n panel.border = element_rect(linewidth = 0.5)\n )","metadata":{"execution":{"iopub.status.busy":"2024-06-29T05:09:42.959656Z","iopub.execute_input":"2024-06-29T05:09:42.961427Z","iopub.status.idle":"2024-06-29T05:09:43.365822Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(25, 15)\n\nggplot(bb_summary, aes(x = P_act_int)) +\n geom_bar(stat = \"count\", color = \"black\", fill = \"steelblue\") +\n geom_text(stat = \"count\", aes(label = after_stat(count)), vjust = -0.5, size = 6) +\n scale_y_continuous(expand = expansion(mult = c(0, 0.3))) +\n facet_wrap(~ protein_name + BB, scales = \"free_y\") +\n ggtitle(\"Numbers of BB at different P(active) intervals across all the 3 proteins and 3 positions\") +\n xlab(\"P(active) interval\") +\n theme_bw(base_size = 22) +\n theme(legend.position = \"bottom\",\n panel.grid.major = element_line(linewidth = 0.2),\n panel.grid.minor = element_line(linewidth = 0.2),\n panel.border = element_rect(linewidth = 0.5),\n strip.background = element_blank()\n )","metadata":{"execution":{"iopub.status.busy":"2024-06-29T05:10:08.841472Z","iopub.execute_input":"2024-06-29T05:10:08.843279Z","iopub.status.idle":"2024-06-29T05:10:10.779438Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Proportion of BB in each P(active) interval across different positions and proteins**","metadata":{}},{"cell_type":"code","source":"bb_summary[, n_gr := .N, by = c(\"protein_name\", \"BB\")]\n\nfreq_dt <- unique(bb_summary[, .(n = .N, freq = round(.N/n_gr, 3)),\n by = c(\"protein_name\", \"BB\", \"P_act_int\")])\ndcast(freq_dt[, !\"n\"], ... ~ protein_name, value.var = \"freq\", fill = 0)","metadata":{"execution":{"iopub.status.busy":"2024-06-29T05:10:32.371729Z","iopub.execute_input":"2024-06-29T05:10:32.373485Z","iopub.status.idle":"2024-06-29T05:10:32.420496Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**Median P(active) value of all BBs across different positions and proteins**","metadata":{}},{"cell_type":"code","source":"dcast(bb_summary[, .(median_p_active = median(p_active)),\n by = c(\"protein_name\", \"BB\")],\n ... ~ protein_name, value.var = \"median_p_active\", fill = 0\n)","metadata":{"execution":{"iopub.status.busy":"2024-06-29T05:11:45.403030Z","iopub.execute_input":"2024-06-29T05:11:45.405059Z","iopub.status.idle":"2024-06-29T05:11:45.438384Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Do the smiles shared across BB2 and BB3 have different P(active) score?","metadata":{}},{"cell_type":"code","source":"bb_summary[, c(\"P_act_int\", \"n_gr\") := NULL]\n\nshared_smiles <- intersect(\n bb_summary[BB == \"buildingblock2_smiles\", SMILES],\n bb_summary[BB == \"buildingblock3_smiles\", SMILES]\n)\nlength(shared_smiles)\n\nshared_smiles_dt <- bb_summary[SMILES %in% shared_smiles]","metadata":{"execution":{"iopub.status.busy":"2024-06-29T05:12:16.128608Z","iopub.execute_input":"2024-06-29T05:12:16.130273Z","iopub.status.idle":"2024-06-29T05:12:16.182065Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(25, 8)\n\nggplot(shared_smiles_dt, aes(x = p_active, fill = BB)) +\n geom_histogram(bins = 30, position = \"dodge\") +\n labs(title = \"Histograms of P(active) Scores for BB2 and BB3\") +\n scale_fill_manual(values = c(\"buildingblock2_smiles\" = \"steelblue\",\n \"buildingblock3_smiles\" = \"#f0cb67\")) +\n facet_wrap(~ protein_name, scales = \"free\") +\n theme_bw(base_size = 22) +\n theme(legend.position = \"bottom\",\n panel.grid.major = element_line(linewidth = 0.2),\n panel.grid.minor = element_line(linewidth = 0.2),\n panel.border = element_rect(linewidth = 0.5),\n strip.background = element_blank()\n )","metadata":{"execution":{"iopub.status.busy":"2024-06-29T05:12:26.037862Z","iopub.execute_input":"2024-06-29T05:12:26.039589Z","iopub.status.idle":"2024-06-29T05:12:26.843244Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"shared_smiles_dt <- dcast(\n shared_smiles_dt, protein_name + SMILES ~ BB,\n value.var = \"p_active\", fill = 0\n )[order(SMILES)]\n\ntst <- shared_smiles_dt[\n , .(p_val = wilcox.test(buildingblock2_smiles, buildingblock3_smiles,\n paired = TRUE)$p.value),\n by = protein_name\n]\n\ntst[, p_val := ifelse(p_val < 0.001, \"<0.001\", round(p_val, 3))]\ntst","metadata":{"execution":{"iopub.status.busy":"2024-06-29T05:12:35.366044Z","iopub.execute_input":"2024-06-29T05:12:35.367754Z","iopub.status.idle":"2024-06-29T05:12:35.417054Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"shared_smiles_dt[\n , p_active_fold_diff := buildingblock3_smiles/buildingblock2_smiles\n]\nshared_smiles_dt[p_active_fold_diff > 2, .N] / shared_smiles_dt[, .N]","metadata":{"execution":{"iopub.status.busy":"2024-06-29T05:12:57.475919Z","iopub.execute_input":"2024-06-29T05:12:57.477625Z","iopub.status.idle":"2024-06-29T05:12:57.498908Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"rm(breaks, freq_dt, max_active, shared_smiles, shared_smiles_dt) ; g <- gc()","metadata":{"execution":{"iopub.status.busy":"2024-06-29T05:13:06.366816Z","iopub.execute_input":"2024-06-29T05:13:06.368532Z","iopub.status.idle":"2024-06-29T05:13:06.591288Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### What characterizes a high P(active) building block?","metadata":{}},{"cell_type":"markdown","source":"First I extract from the train data set all molecules with binds = 1 (molecules that bind to at least one protein)","metadata":{}},{"cell_type":"code","source":"bind1 <- dbGetQuery(con, paste0(\n\"SELECT * FROM train\nWHERE binds = 1\"\n))\nsetDT(bind1)","metadata":{"execution":{"iopub.status.busy":"2024-06-29T05:13:14.350249Z","iopub.execute_input":"2024-06-29T05:13:14.351948Z","iopub.status.idle":"2024-06-29T05:14:38.194918Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Btw, are there molecules, where building blocks 2 and are the same?","metadata":{}},{"cell_type":"code","source":"bind1[buildingblock2_smiles == buildingblock3_smiles, .N]","metadata":{"execution":{"iopub.status.busy":"2024-04-21T08:58:02.839081Z","iopub.execute_input":"2024-04-21T08:58:02.840745Z","iopub.status.idle":"2024-04-21T08:58:02.871568Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Building blocks with a higher number of different building blocks in molecules that bind to proteins may indicate a propensity to bind. So I subset BBs that are present in at least 1 binding molecule. For each of these BBs, I calculate the number of other BBs that co-occur with a given BB. I use this to see if there is a relationship between the number of co-occurring other BBs and the P(active) score.","metadata":{}},{"cell_type":"code","source":"bind1 <- melt(bind1, measure.vars = bbs_cols,\n value.name = \"SMILES\", variable.name = \"BB\",\n variable.factor = FALSE)\n\nbind1 <- bind1[, .(molecule_smiles, protein_name, BB, SMILES)]\nactive_bbs <- bb_summary[num_active > 0, .(protein_name, BB, SMILES, p_active)]","metadata":{"execution":{"iopub.status.busy":"2024-04-21T08:58:02.874429Z","iopub.execute_input":"2024-04-21T08:58:02.875968Z","iopub.status.idle":"2024-04-21T08:58:08.992073Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tictoc::tic()\ntop_bbs_mols <- rbindlist(\n lapply(1:active_bbs[, .N], function(r) {\n \n mols <- bind1[protein_name == active_bbs[r, protein_name] &\n BB == active_bbs[r, BB] & \n SMILES == active_bbs[r, SMILES], unique(molecule_smiles)]\n \n res <- bind1[molecule_smiles %chin% mols &\n protein_name == active_bbs[r, protein_name]]\n \n res <- res[, .(n_unique_smiles = uniqueN(SMILES)),\n by = c(\"protein_name\", \"BB\")]\n \n res <- dcast(res, ... ~ BB, value.var = \"n_unique_smiles\")\n \n res[, `:=` (target_smiles = active_bbs[r, SMILES],\n target_BB = active_bbs[r, BB],\n p_active = active_bbs[r, p_active]\n )\n ]\n res\n })\n)\ntictoc::toc()","metadata":{"execution":{"iopub.status.busy":"2024-04-21T08:58:08.999160Z","iopub.execute_input":"2024-04-21T08:58:09.002692Z","iopub.status.idle":"2024-04-21T09:02:33.330631Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"top_bbs_mols <- melt(top_bbs_mols,\n measure.vars = bbs_cols,\n value.name = \"n_distinct_bbs\",\n variable.name = \"other_BB\",\n variable.factor = FALSE)\nfreqs <- bb_smiles[, .(n_BB_smiles = .N), by = \"BB\"]\n\ntop_bbs_mols <- top_bbs_mols[freqs, on = c(\"other_BB\" = \"BB\")]\n\n\ntop_bbs_mols[\n , prop_distinct_bbs := ifelse(\n other_BB != target_BB, n_distinct_bbs / n_BB_smiles, n_distinct_bbs\n )\n]","metadata":{"execution":{"iopub.status.busy":"2024-04-21T09:02:33.335186Z","iopub.execute_input":"2024-04-21T09:02:33.337038Z","iopub.status.idle":"2024-04-21T09:02:33.369296Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig(25, 20)\n\nggplot(top_bbs_mols[target_BB != other_BB],\n aes(x = p_active, y = prop_distinct_bbs)) +\n geom_point(shape = 21, size = 0.7) +\n theme_bw(base_size = 22) +\n theme(panel.grid.major = element_blank(),\n panel.grid.minor = element_blank()) +\n facet_grid(target_BB ~ other_BB, drop = TRUE) +\n theme(strip.placement = \"outside\",\n strip.background = element_blank())","metadata":{"execution":{"iopub.status.busy":"2024-04-21T09:02:33.372341Z","iopub.execute_input":"2024-04-21T09:02:33.374032Z","iopub.status.idle":"2024-04-21T09:02:35.329352Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Adding log scale for x axis","metadata":{}},{"cell_type":"code","source":"ggplot(top_bbs_mols[target_BB != other_BB],\n aes(x = p_active, y = prop_distinct_bbs)) +\n geom_point(shape = 21, size = 0.7) +\n theme_bw(base_size = 22) +\n theme(panel.grid.major = element_blank(),\n panel.grid.minor = element_blank()) +\n facet_grid(target_BB ~ other_BB, drop = TRUE) +\n ggtitle(\"P(active) vs. Proportion of Distinct BBs across different positions\") +\n theme(strip.placement = \"outside\",\n strip.background = element_blank()) +\n scale_x_continuous(trans = \"log10\")","metadata":{"execution":{"iopub.status.busy":"2024-04-21T09:02:35.332846Z","iopub.execute_input":"2024-04-21T09:02:35.334688Z","iopub.status.idle":"2024-04-21T09:02:37.241335Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"From the article: P(active) tells us how successful a building block is when it is placed in a certain position, but tells us nothing about the behavior at other positions that may lead to the success (or lack thereof) at one position. Hypothetically, a building block could have a high P(active) value but only form active compounds with a very limited selection of partners in one of the other positions. On the contrary, we see that building blocks which are successful in one position are compatible with a broader diversity of building blocks in all other positions. \n\nI think we see a similar picture here.","metadata":{}}]}