{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Adding Compound Info from Optional Files\n\nThis [post](https://www.kaggle.com/competitions/open-problems-single-cell-perturbations/discussion/440680) linked an s3 bucket with some\noptional files: `s3://saturn-kaggle-datasets/open-problems-single-cell-perturbations-optional/`\n\nThis notebooks works through reading the `compoundinfo_beta.txt` file from there and aligning with our `sm_names` field using either string matching or `SMILES` based matching. We then explore the mechanism of action and target information in the dimensionally reduced UMAP embeddings of the gene expression space.","metadata":{}},{"cell_type":"code","source":"import pandas as pd\nfrom pathlib import Path","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-09-24T03:26:34.383130Z","iopub.execute_input":"2023-09-24T03:26:34.384082Z","iopub.status.idle":"2023-09-24T03:26:34.390136Z","shell.execute_reply.started":"2023-09-24T03:26:34.384039Z","shell.execute_reply":"2023-09-24T03:26:34.388159Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"I've loaded the file up to Kaggle as a dataset and access it in this notebook as additional mounted data with `+ Add Data`.","metadata":{}},{"cell_type":"code","source":"data_dir1 = Path(\"/kaggle/input/open-problems-single-cell-perturbations\")\ndata_dir2 = Path(\"/kaggle/input/supplemental-compound-info-for-op2\")\nchem_info = pd.read_csv(data_dir2 / \"compoundinfo_beta.txt\", sep='\\t')\nde_train = pd.read_parquet(data_dir1 / \"de_train.parquet\")","metadata":{"execution":{"iopub.status.busy":"2023-09-24T03:16:25.702258Z","iopub.execute_input":"2023-09-24T03:16:25.702665Z","iopub.status.idle":"2023-09-24T03:16:27.095054Z","shell.execute_reply.started":"2023-09-24T03:16:25.702634Z","shell.execute_reply":"2023-09-24T03:16:27.094072Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Next we use our `sm_names` and `smiles` fields to look for matches in the compoundinfo txt (TSV) file.","metadata":{}},{"cell_type":"code","source":"sm_names = list(de_train['sm_name'].unique())\nsmiles = list(de_train['SMILES'].unique())","metadata":{"execution":{"iopub.status.busy":"2023-09-24T03:26:39.021677Z","iopub.execute_input":"2023-09-24T03:26:39.022085Z","iopub.status.idle":"2023-09-24T03:26:39.028704Z","shell.execute_reply.started":"2023-09-24T03:26:39.022054Z","shell.execute_reply":"2023-09-24T03:26:39.027512Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"smile_set = set(chem_info['canonical_smiles'])\nfound_smile = [smile for smile in smiles if smile in smile_set]\nprint(f\"Matched: {len(found_smile)} out of {len(smiles)} molecules by smile representation.\")","metadata":{"execution":{"iopub.status.busy":"2023-09-24T03:26:39.620954Z","iopub.execute_input":"2023-09-24T03:26:39.621801Z","iopub.status.idle":"2023-09-24T03:26:39.637937Z","shell.execute_reply.started":"2023-09-24T03:26:39.621762Z","shell.execute_reply":"2023-09-24T03:26:39.636646Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"found_smname = [sm_name.upper() for sm_name in sm_names\n                if sm_name.upper() in set(chem_info['cmap_name'].str.upper())]\nprint(f\"Matched: {len(found_smname)} out of {len(sm_names)} molecules by name.\")","metadata":{"execution":{"iopub.status.busy":"2023-09-24T03:26:43.188366Z","iopub.execute_input":"2023-09-24T03:26:43.188782Z","iopub.status.idle":"2023-09-24T03:26:46.279912Z","shell.execute_reply.started":"2023-09-24T03:26:43.188751Z","shell.execute_reply":"2023-09-24T03:26:46.278707Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We subset both our `de_train` and `chem_info` dataframes to contain only rows we know will align. This is to simplify further exploration. For the best results in this competition, you probably want to manually annotate unmatched compounds and/or look for other databases where smiles or LINCS ID can be used to locate targets and mechanism of action information.","metadata":{}},{"cell_type":"code","source":"match_df = de_train[de_train['SMILES'].isin(found_smile) |\n                    de_train['sm_name'].str.upper().isin(found_smname)]\nmatch_df = match_df.copy(deep=True)\nmatch_df.head(10)","metadata":{"execution":{"iopub.status.busy":"2023-09-24T03:26:46.282055Z","iopub.execute_input":"2023-09-24T03:26:46.282728Z","iopub.status.idle":"2023-09-24T03:26:46.494343Z","shell.execute_reply.started":"2023-09-24T03:26:46.282683Z","shell.execute_reply":"2023-09-24T03:26:46.493306Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Just under 72% of the rows containing compounds can be matched.","metadata":{}},{"cell_type":"code","source":"len(match_df) / len(de_train)","metadata":{"execution":{"iopub.status.busy":"2023-09-24T03:26:49.631925Z","iopub.execute_input":"2023-09-24T03:26:49.632395Z","iopub.status.idle":"2023-09-24T03:26:49.640238Z","shell.execute_reply.started":"2023-09-24T03:26:49.632359Z","shell.execute_reply":"2023-09-24T03:26:49.638980Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"chem_df = chem_info[chem_info['cmap_name'].str.upper().isin(found_smname) |\n                    chem_info['canonical_smiles'].isin(found_smile)]\nchem_df = chem_df.copy(deep=True)\nchem_df.head(10)","metadata":{"execution":{"iopub.status.busy":"2023-09-24T03:26:50.482277Z","iopub.execute_input":"2023-09-24T03:26:50.482716Z","iopub.status.idle":"2023-09-24T03:26:50.522627Z","shell.execute_reply.started":"2023-09-24T03:26:50.482682Z","shell.execute_reply":"2023-09-24T03:26:50.521336Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Better news: less than 5% of our aligned compound info is missing mechanism of action information.","metadata":{}},{"cell_type":"code","source":"sum(chem_df['moa'].isna()) / len(chem_df['moa'])","metadata":{"execution":{"iopub.status.busy":"2023-09-24T03:26:53.297687Z","iopub.execute_input":"2023-09-24T03:26:53.298122Z","iopub.status.idle":"2023-09-24T03:26:53.306518Z","shell.execute_reply.started":"2023-09-24T03:26:53.298089Z","shell.execute_reply":"2023-09-24T03:26:53.305155Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Joining Dataframes\n\nWe join both directions, for different purposes.\n\nFirst we create a field with a uniform uppercase representation of the compound name.","metadata":{}},{"cell_type":"code","source":"chem_df['synth_id'] = chem_df['cmap_name'].str.upper()\nmatch_df['synth_id'] = match_df['sm_name'].str.upper()\naligned = match_df.merge(chem_df[['synth_id', 'cmap_name']], how=\"left\", on='synth_id')","metadata":{"execution":{"iopub.status.busy":"2023-09-24T03:26:55.547755Z","iopub.execute_input":"2023-09-24T03:26:55.548194Z","iopub.status.idle":"2023-09-24T03:26:55.940231Z","shell.execute_reply.started":"2023-09-24T03:26:55.548160Z","shell.execute_reply":"2023-09-24T03:26:55.939005Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We note that 26 remain unmatched. Since we limited this to either matching name or matching SMILES information, this means we need to align using SMILES.","metadata":{}},{"cell_type":"code","source":"unmatched = aligned[aligned['cmap_name'].isna()]\nlen(unmatched)","metadata":{"execution":{"iopub.status.busy":"2023-09-24T03:26:58.389396Z","iopub.execute_input":"2023-09-24T03:26:58.389814Z","iopub.status.idle":"2023-09-24T03:26:58.403236Z","shell.execute_reply.started":"2023-09-24T03:26:58.389784Z","shell.execute_reply":"2023-09-24T03:26:58.401734Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now we verify that we can match everything we linked previously using SMILES in the case that we can't link it by ID. We inspect the resulting matched compounds (can verify this via google searches, online DBs, that these are the same compounds).","metadata":{}},{"cell_type":"code","source":"test = unmatched.merge(chem_df, left_on='SMILES', right_on='canonical_smiles', how='left')\nprint(f\"Remaining unmatched: {sum(test['synth_id_y'].isna())}.\")\ntest[['sm_name', 'cmap_name_y']].drop_duplicates()","metadata":{"execution":{"iopub.status.busy":"2023-09-24T03:27:01.082535Z","iopub.execute_input":"2023-09-24T03:27:01.082997Z","iopub.status.idle":"2023-09-24T03:27:01.142153Z","shell.execute_reply.started":"2023-09-24T03:27:01.082964Z","shell.execute_reply":"2023-09-24T03:27:01.141014Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"compounds = chem_df.merge(unmatched[['SMILES', 'sm_name']], right_on='SMILES', left_on='canonical_smiles', how='left')\ncompounds = compounds.merge(match_df[['sm_name', 'synth_id']], how='left', on='synth_id')\ncompounds.head()","metadata":{"execution":{"iopub.status.busy":"2023-09-24T03:27:03.992964Z","iopub.execute_input":"2023-09-24T03:27:03.993601Z","iopub.status.idle":"2023-09-24T03:27:04.018176Z","shell.execute_reply.started":"2023-09-24T03:27:03.993570Z","shell.execute_reply":"2023-09-24T03:27:04.017035Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Still a good ratio of not MOA NA values.","metadata":{}},{"cell_type":"code","source":"sum(compounds['moa'].isna()) / len(compounds)","metadata":{"execution":{"iopub.status.busy":"2023-09-24T03:27:07.876587Z","iopub.execute_input":"2023-09-24T03:27:07.877031Z","iopub.status.idle":"2023-09-24T03:27:07.884788Z","shell.execute_reply.started":"2023-09-24T03:27:07.876995Z","shell.execute_reply":"2023-09-24T03:27:07.883621Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now let's clean up the compounds DataFrame.","metadata":{}},{"cell_type":"code","source":"compounds.columns","metadata":{"execution":{"iopub.status.busy":"2023-09-24T03:27:10.470345Z","iopub.execute_input":"2023-09-24T03:27:10.470760Z","iopub.status.idle":"2023-09-24T03:27:10.478730Z","shell.execute_reply.started":"2023-09-24T03:27:10.470727Z","shell.execute_reply":"2023-09-24T03:27:10.477497Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"compounds['sm_name'] = compounds['sm_name_x'].fillna(compounds['sm_name_y'])\ncompounds.drop(columns=['sm_name_x', 'sm_name_y'], inplace=True)\ncompounds.head(8)","metadata":{"execution":{"iopub.status.busy":"2023-09-24T03:27:11.172900Z","iopub.execute_input":"2023-09-24T03:27:11.173643Z","iopub.status.idle":"2023-09-24T03:27:11.196541Z","shell.execute_reply.started":"2023-09-24T03:27:11.173596Z","shell.execute_reply":"2023-09-24T03:27:11.195739Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"compounds.to_csv('/kaggle/working/compounds.tsv', sep='\\t')","metadata":{"execution":{"iopub.status.busy":"2023-09-24T03:27:14.006857Z","iopub.execute_input":"2023-09-24T03:27:14.007260Z","iopub.status.idle":"2023-09-24T03:27:14.036016Z","shell.execute_reply.started":"2023-09-24T03:27:14.007230Z","shell.execute_reply":"2023-09-24T03:27:14.035041Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now let's look at a breakdown of mechanism of action to look for a few classes we might be able to visualize.","metadata":{}},{"cell_type":"code","source":"from matplotlib import pyplot as plt\n\nmoa_agg = compounds['moa'].value_counts()\nprint(moa_agg.iloc[:10])\nmoa_agg.hist()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-09-24T03:27:15.300715Z","iopub.execute_input":"2023-09-24T03:27:15.301158Z","iopub.status.idle":"2023-09-24T03:27:15.554832Z","shell.execute_reply.started":"2023-09-24T03:27:15.301123Z","shell.execute_reply":"2023-09-24T03:27:15.553725Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%capture\n# ^ capture is to make it shut up about umap numba warnings, I don't\n#   actually use the captured value.\nfrom umap.umap_ import UMAP\nimport plotly.express as px\nfrom sklearn.decomposition import TruncatedSVD","metadata":{"execution":{"iopub.status.busy":"2023-09-24T03:29:26.020494Z","iopub.execute_input":"2023-09-24T03:29:26.020920Z","iopub.status.idle":"2023-09-24T03:29:27.077765Z","shell.execute_reply.started":"2023-09-24T03:29:26.020888Z","shell.execute_reply":"2023-09-24T03:29:27.076497Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We've got a bit of a tricky issue in the join, because compounds has multiple rows, mostly to handle multiple targets, but also several compounds have more than one mechanism of action.\n\nThe below oesn't work because it's still not enough to remove duplicate sm_name or cmap_name values. So\nwe go back to naive geneX space and join on info later.","metadata":{}},{"cell_type":"code","source":"many_targets = compounds.groupby('sm_name')['target'].apply(lambda x: list(set(x)))\nmany_moa = compounds.groupby('sm_name')['target'].apply(lambda x: list(set(x)))\nsimple_compounds = compounds.drop(columns=['target', 'compound_aliases', 'moa']).drop_duplicates()\n# simple_compounds['targets'] = many_targets.values\n# simple_compounds['moas'] = many_moa.values\nsimple_compounds.head(1)","metadata":{"execution":{"iopub.status.busy":"2023-09-24T03:44:19.934800Z","iopub.execute_input":"2023-09-24T03:44:19.935245Z","iopub.status.idle":"2023-09-24T03:44:19.965418Z","shell.execute_reply.started":"2023-09-24T03:44:19.935208Z","shell.execute_reply":"2023-09-24T03:44:19.964175Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Visualize Compound Info in Gene Expression Space\n\nEven though we only have limited info to join, we need original geneX space for UMAP projection.","metadata":{}},{"cell_type":"code","source":"gX = de_train[de_train.columns[5:]].values.copy()\ntsvd = TruncatedSVD(n_components=35, n_iter=7, random_state=10)\ngX_dr = tsvd.fit_transform(gX)\nump = UMAP(n_components=3)\ngX_proj = ump.fit_transform(gX_dr)","metadata":{"execution":{"iopub.status.busy":"2023-09-24T03:46:32.832398Z","iopub.execute_input":"2023-09-24T03:46:32.833491Z","iopub.status.idle":"2023-09-24T03:46:42.946428Z","shell.execute_reply.started":"2023-09-24T03:46:32.833450Z","shell.execute_reply":"2023-09-24T03:46:42.945457Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We need to put UMAP X,Y,Z projections in a dataframe subset where we'll join chemical information and add visualization fields.","metadata":{}},{"cell_type":"code","source":"viz_df = de_train.drop(columns=de_train.columns[5:]).copy(deep=True)\nviz_df['umapX'] = gX_proj[:, 0]\nviz_df['umapY'] = gX_proj[:, 1]\nviz_df['umapZ'] = gX_proj[:, 2]\nviz_df.head(5)","metadata":{"execution":{"iopub.status.busy":"2023-09-24T03:53:32.207382Z","iopub.execute_input":"2023-09-24T03:53:32.207937Z","iopub.status.idle":"2023-09-24T03:53:32.228217Z","shell.execute_reply.started":"2023-09-24T03:53:32.207893Z","shell.execute_reply":"2023-09-24T03:53:32.227319Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"viz_df['synth_id'] = viz_df['sm_name'].str.upper()\nviz_df = viz_df.merge(compounds, on='synth_id', how='left')","metadata":{"execution":{"iopub.status.busy":"2023-09-24T03:49:12.317246Z","iopub.execute_input":"2023-09-24T03:49:12.317713Z","iopub.status.idle":"2023-09-24T03:49:12.333057Z","shell.execute_reply.started":"2023-09-24T03:49:12.317678Z","shell.execute_reply":"2023-09-24T03:49:12.331848Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"few_classes = [\n    'VEGFR inhibitor',\n    'PDGFR inhibitor',\n    'CDK inhibitor',\n    'KIT inhibitor',\n    'FLT3 inhibitor',\n]\nviz_df['narrow_class'] = viz_df['moa'].apply(lambda x: x if x in few_classes else 'Other')\nviz_df['size_hint'] = viz_df['moa'].apply(lambda x: 8 if x in few_classes else 2)","metadata":{"execution":{"iopub.status.busy":"2023-09-24T03:51:14.823047Z","iopub.execute_input":"2023-09-24T03:51:14.823511Z","iopub.status.idle":"2023-09-24T03:51:14.840130Z","shell.execute_reply.started":"2023-09-24T03:51:14.823476Z","shell.execute_reply":"2023-09-24T03:51:14.838830Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"px.scatter_3d(\n    viz_df,\n    x='umapX',\n    y='umapY',\n    z='umapZ',\n    color='narrow_class',\n    size='size_hint',\n    hover_data=['sm_name_y', 'moa', 'target']\n)","metadata":{"execution":{"iopub.status.busy":"2023-09-24T03:51:19.393330Z","iopub.execute_input":"2023-09-24T03:51:19.393752Z","iopub.status.idle":"2023-09-24T03:51:19.587739Z","shell.execute_reply.started":"2023-09-24T03:51:19.393716Z","shell.execute_reply":"2023-09-24T03:51:19.586763Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Is it possible cell type is a confounding variable and we just need to look at these MoAs on one cell type?","metadata":{}},{"cell_type":"code","source":"de_train['cell_type'].unique()","metadata":{"execution":{"iopub.status.busy":"2023-09-24T03:54:12.930597Z","iopub.execute_input":"2023-09-24T03:54:12.930977Z","iopub.status.idle":"2023-09-24T03:54:12.939010Z","shell.execute_reply.started":"2023-09-24T03:54:12.930946Z","shell.execute_reply":"2023-09-24T03:54:12.937602Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We repeat our prep steps and UMAP embeddings for CD4+ cells only.","metadata":{}},{"cell_type":"code","source":"cd4 = de_train[de_train['cell_type'] == 'T cells CD4+']\ngX = cd4[cd4.columns[5:]].values.copy()\ntsvd = TruncatedSVD(n_components=35, n_iter=7, random_state=10)\ngX_dr = tsvd.fit_transform(gX)\nump = UMAP(n_components=3)\ngX_proj = ump.fit_transform(gX_dr)\ngX_proj.shape","metadata":{"execution":{"iopub.status.busy":"2023-09-24T03:56:41.133812Z","iopub.execute_input":"2023-09-24T03:56:41.134628Z","iopub.status.idle":"2023-09-24T03:56:45.941083Z","shell.execute_reply.started":"2023-09-24T03:56:41.134584Z","shell.execute_reply":"2023-09-24T03:56:45.939810Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We then join all visualization information into one data frame, with gene expression removed.","metadata":{}},{"cell_type":"code","source":"viz_df = cd4.drop(columns=cd4.columns[5:]).copy(deep=True)\nviz_df['umapX'] = gX_proj[:, 0]\nviz_df['umapY'] = gX_proj[:, 1]\nviz_df['umapZ'] = gX_proj[:, 2]\nviz_df['synth_id'] = viz_df['sm_name'].str.upper()\nviz_df = viz_df.merge(compounds, on='synth_id', how='left')\nviz_df['narrow_class'] = viz_df['moa'].apply(lambda x: x if x in few_classes else 'Other')\nviz_df['size_hint'] = viz_df['moa'].apply(lambda x: 8 if x in few_classes else 2)","metadata":{"execution":{"iopub.status.busy":"2023-09-24T03:56:48.732338Z","iopub.execute_input":"2023-09-24T03:56:48.732790Z","iopub.status.idle":"2023-09-24T03:56:48.754635Z","shell.execute_reply.started":"2023-09-24T03:56:48.732757Z","shell.execute_reply":"2023-09-24T03:56:48.753451Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"px.scatter_3d(\n    viz_df,\n    x='umapX',\n    y='umapY',\n    z='umapZ',\n    color='narrow_class',\n    size='size_hint',\n    hover_data=['sm_name_y', 'moa', 'target']\n)","metadata":{"execution":{"iopub.status.busy":"2023-09-24T03:56:53.688871Z","iopub.execute_input":"2023-09-24T03:56:53.690072Z","iopub.status.idle":"2023-09-24T03:56:53.818244Z","shell.execute_reply.started":"2023-09-24T03:56:53.690030Z","shell.execute_reply":"2023-09-24T03:56:53.816969Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"However even without the change, it's not clear that our most common inhibitor classes seem to be in similar places in gene expression space. We have more work to do to figure out how to characterize chemicals in a way that looks like it aligns with gene expression outcomes in perturbation experiments.\n\nThis notebook is the work of Vendekagon Labs and licensed under Apache 2.0, as all public Kaggle notebooks are. If you found this helpful, you may also get something out of:\n\n- [Pathway Enrichment Analysis](https://www.kaggle.com/code/vendekagonlabs/op2-pathway-enrichment)\n- [Single Cell Analysis with Scanpy](https://www.kaggle.com/code/vendekagonlabs/op2-single-cell-eda-10x-multiome)\n- [Chemical Feature Exploration in R](https://www.kaggle.com/code/vendekagonlabs/come-on-chemicals-r-version)","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}