{"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":"# Pathway / Gene Set Enrichment Analysis\n\nThis notebook uses GSEApy to explore using pathway enrichment on our gene\nexpression data. Pathway/Gene set enrichment associates genes with known\npathways or related genes that share a common function or metabolic path, etc.\nThis would let us take a set of gene symbols that are up or down-regulated in\nreaction to treatment with a molecule, and see which pathways are impacted.\n\nFor example, it would be simpler to see that Apoptosis pathways are upregulated and\nproliferation pathways are downregulated using prior knowledge encoded in gene sets,\nrather than trying to piece this together ourselves by examining every gene's function\nin isolation.","metadata":{}},{"cell_type":"code","source":"%%capture\n%pip install scanpy gseapy leidenalg python-igraph","metadata":{"execution":{"iopub.status.busy":"2023-09-24T15:30:28.305071Z","iopub.execute_input":"2023-09-24T15:30:28.305474Z","iopub.status.idle":"2023-09-24T15:30:40.751012Z","shell.execute_reply.started":"2023-09-24T15:30:28.305441Z","shell.execute_reply":"2023-09-24T15:30:40.749599Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import pandas as pd\nimport scanpy as sc\nimport gseapy as gp\nfrom matplotlib import pyplot as plt\nfrom pathlib import Path\nfrom time import sleep\n\ndata_dir = Path('/kaggle/input/open-problems-single-cell-perturbations')\nde_train = pd.read_parquet(data_dir / 'de_train.parquet')\nde_train.head(10)","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-09-24T15:30:40.753948Z","iopub.execute_input":"2023-09-24T15:30:40.754274Z","iopub.status.idle":"2023-09-24T15:30:42.022681Z","shell.execute_reply.started":"2023-09-24T15:30:40.754246Z","shell.execute_reply":"2023-09-24T15:30:42.021630Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Simple Walkthrough\n\nFor a simple walkthrough of basic enrichment analysis, we isolate to one cell type and two molecules with different functions.","metadata":{}},{"cell_type":"code","source":"cell_type = 'T cells CD4+'\ngenes = de_train.columns[5:]\n# Idelalisib - PI3K inhibitor (Leukemia & Lymphoma): prevents proliferation\n#.             & induces apoptosis\n# Alogliptin - DPP-4 inhibitor (Diabetes): lowers blood sugar\nmolecules = ['Idelalisib', 'Alogliptin']\n\nsub_df = de_train[(de_train['sm_name'].isin(molecules)) &\n                  (de_train['cell_type'] == cell_type)]\nsub_df","metadata":{"execution":{"iopub.status.busy":"2023-09-24T15:30:42.023830Z","iopub.execute_input":"2023-09-24T15:30:42.024123Z","iopub.status.idle":"2023-09-24T15:30:42.052535Z","shell.execute_reply.started":"2023-09-24T15:30:42.024097Z","shell.execute_reply":"2023-09-24T15:30:42.051640Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Heuristic based gene selection\n\nThe most basic enrichment analysis lets you supply a set of genes using your own criteria, then informs you about whether or not the genes you provide show significant overlap with a known pathway. In this case, we just select the most highly upregulated and highly downregulated genes (lowest and highest expression change).","metadata":{}},{"cell_type":"code","source":"idelalisib_expr = sub_df[genes].loc[11]\nq0, q1 = idelalisib_expr.quantile(0.01), idelalisib_expr.quantile(0.99)\nprint(f\"cutoffs: {q0}, {q1}\")\n\nfig, ax = plt.subplots(1, 1)\nidelalisib_expr.hist(bins=50, ax=ax)\nax.set_ylim(0, 5000)\nx1, y1 = [q0, q0], [0, 5000]\nx2, y2 = [q1, q1], [0, 5000]\nax.plot(x1, y1, x2, y2, marker = '.', linestyle='--')\nax.set_title('0.01 and 0.99 logfold expression cutoffs.')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-09-24T15:30:42.054671Z","iopub.execute_input":"2023-09-24T15:30:42.054947Z","iopub.status.idle":"2023-09-24T15:30:42.418723Z","shell.execute_reply.started":"2023-09-24T15:30:42.054923Z","shell.execute_reply":"2023-09-24T15:30:42.417658Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"idelalisib_up = sub_df[genes].loc[11].sort_values(ascending=False)[:500]\nidelalisib_down = sub_df[genes].loc[11].sort_values(ascending=True)[:500]\nalogliptin_up = sub_df[genes].loc[464].sort_values(ascending=False)[:500]\nalogliptin_down = sub_df[genes].loc[464].sort_values(ascending=True)[:500]","metadata":{"execution":{"iopub.status.busy":"2023-09-24T15:30:42.420738Z","iopub.execute_input":"2023-09-24T15:30:42.421063Z","iopub.status.idle":"2023-09-24T15:30:42.457816Z","shell.execute_reply.started":"2023-09-24T15:30:42.421034Z","shell.execute_reply":"2023-09-24T15:30:42.457008Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def enrichr_query(gene_list):\n    return gp.enrichr(gene_list=gene_list,\n                      gene_sets=['KEGG_2021_Human', 'GO_Biological_Process_2021'],\n                      organism='human',\n                      outdir=None)\n\n\ndef enrichr_results(gene_list, retry_pace=8):\n    while True:\n        try:\n            result = enrichr_query(gene_list)\n            break\n        except:\n            print(f\"Server failed to respond, retrying in {retry_pace} seconds.\")\n            sleep(retry_pace)\n    return result","metadata":{"execution":{"iopub.status.busy":"2023-09-24T15:30:42.458797Z","iopub.execute_input":"2023-09-24T15:30:42.459657Z","iopub.status.idle":"2023-09-24T15:30:42.465942Z","shell.execute_reply.started":"2023-09-24T15:30:42.459622Z","shell.execute_reply":"2023-09-24T15:30:42.464961Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"gene_list = list(idelalisib_up[idelalisib_up > q1].index)\nid_up = enrichr_results(gene_list)\nid_up.results[id_up.results['Adjusted P-value'] > 0.05]","metadata":{"execution":{"iopub.status.busy":"2023-09-24T15:30:42.467289Z","iopub.execute_input":"2023-09-24T15:30:42.467633Z","iopub.status.idle":"2023-09-24T15:30:45.542334Z","shell.execute_reply.started":"2023-09-24T15:30:42.467559Z","shell.execute_reply":"2023-09-24T15:30:45.541274Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"gene_list = list(idelalisib_down[idelalisib_down < q0].index)\nid_down = enrichr_results(gene_list)\nid_down.results[id_down.results['Adjusted P-value'] > 0.05]","metadata":{"execution":{"iopub.status.busy":"2023-09-24T15:30:45.543840Z","iopub.execute_input":"2023-09-24T15:30:45.544169Z","iopub.status.idle":"2023-09-24T15:30:48.413460Z","shell.execute_reply.started":"2023-09-24T15:30:45.544140Z","shell.execute_reply":"2023-09-24T15:30:48.412339Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"gene_list = list(alogliptin_up[alogliptin_up > 2.5].index)\nalg_up = enrichr_results(gene_list)\nalg_up.results.head(10)","metadata":{"execution":{"iopub.status.busy":"2023-09-24T15:30:48.414837Z","iopub.execute_input":"2023-09-24T15:30:48.416069Z","iopub.status.idle":"2023-09-24T15:30:50.614556Z","shell.execute_reply.started":"2023-09-24T15:30:48.416029Z","shell.execute_reply":"2023-09-24T15:30:50.613669Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"gene_list = list(alogliptin_down[alogliptin_down < -2.5].index)\nalg_down = enrichr_results(gene_list)\nalg_down.results.head(10)","metadata":{"execution":{"iopub.status.busy":"2023-09-24T15:30:50.619122Z","iopub.execute_input":"2023-09-24T15:30:50.619913Z","iopub.status.idle":"2023-09-24T15:30:52.788167Z","shell.execute_reply.started":"2023-09-24T15:30:50.619885Z","shell.execute_reply":"2023-09-24T15:30:52.787122Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def dot_plot(enr_obj, cutoff=0.1):\n    ax = gp.dotplot(\n        enr_obj.results,\n        column=\"Adjusted P-value\",\n        x='Combined Score',\n        size=2,\n        top_term=50,\n        figsize=(3,9),\n        title = \"KEGG\",\n        xticklabels_rot=45,\n        show_ring=True,\n        marker='o',\n        cutoff=cutoff,\n    )\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2023-09-24T15:30:52.789740Z","iopub.execute_input":"2023-09-24T15:30:52.790327Z","iopub.status.idle":"2023-09-24T15:30:52.797775Z","shell.execute_reply.started":"2023-09-24T15:30:52.790289Z","shell.execute_reply":"2023-09-24T15:30:52.796523Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dot_plot(id_up, cutoff=0.25)","metadata":{"execution":{"iopub.status.busy":"2023-09-24T15:30:52.799126Z","iopub.execute_input":"2023-09-24T15:30:52.799445Z","iopub.status.idle":"2023-09-24T15:30:53.526918Z","shell.execute_reply.started":"2023-09-24T15:30:52.799420Z","shell.execute_reply":"2023-09-24T15:30:53.525856Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dot_plot(id_down, cutoff=0.0001)","metadata":{"execution":{"iopub.status.busy":"2023-09-24T15:30:53.528406Z","iopub.execute_input":"2023-09-24T15:30:53.528740Z","iopub.status.idle":"2023-09-24T15:30:54.040910Z","shell.execute_reply.started":"2023-09-24T15:30:53.528712Z","shell.execute_reply":"2023-09-24T15:30:54.039659Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dot_plot(alg_up, cutoff=0.062570)","metadata":{"execution":{"iopub.status.busy":"2023-09-24T15:30:54.042994Z","iopub.execute_input":"2023-09-24T15:30:54.043805Z","iopub.status.idle":"2023-09-24T15:30:55.203676Z","shell.execute_reply.started":"2023-09-24T15:30:54.043751Z","shell.execute_reply":"2023-09-24T15:30:55.202822Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dot_plot(alg_down, cutoff=0.05)","metadata":{"execution":{"iopub.status.busy":"2023-09-24T15:30:55.204988Z","iopub.execute_input":"2023-09-24T15:30:55.205546Z","iopub.status.idle":"2023-09-24T15:30:55.925302Z","shell.execute_reply.started":"2023-09-24T15:30:55.205516Z","shell.execute_reply":"2023-09-24T15:30:55.924496Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# General Purpose Enrichment Check per Molecule\n\nYou can replace `random.choice` with a molecule you're investigating. This looks at\npathway enrichment for that particular molecule, using the method defined above.","metadata":{}},{"cell_type":"code","source":"import random\n\nmolecule = random.choice(list(de_train['sm_name']))\nsub_df = de_train[(de_train['sm_name'] == molecule) &\n                  (de_train['cell_type'] == cell_type)]\ngene_expr = sub_df[genes].iloc[0]\nq0, q1 = gene_expr.quantile(0.01), idelalisib_expr.quantile(0.99)\nprint(f\"{molecule} with cutoffs of {q0} and {q1}\")","metadata":{"execution":{"iopub.status.busy":"2023-09-24T15:30:55.926723Z","iopub.execute_input":"2023-09-24T15:30:55.927262Z","iopub.status.idle":"2023-09-24T15:30:55.942124Z","shell.execute_reply.started":"2023-09-24T15:30:55.927233Z","shell.execute_reply":"2023-09-24T15:30:55.941480Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"gene_list = list(gene_expr[gene_expr > q1].index)\nif len(gene_list) > 150:\n    gene_list = gene_list[:100]\nupreg = enrichr_results(gene_list)\nupreg.results.head(10)","metadata":{"execution":{"iopub.status.busy":"2023-09-24T15:30:55.943618Z","iopub.execute_input":"2023-09-24T15:30:55.944397Z","iopub.status.idle":"2023-09-24T15:30:58.606174Z","shell.execute_reply.started":"2023-09-24T15:30:55.944349Z","shell.execute_reply":"2023-09-24T15:30:58.605388Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"try:\n    dot_plot(upreg, cutoff=0.10)\nexcept ValueError:\n    print(\"No pathway match was significant for this gene-set!\")","metadata":{"execution":{"iopub.status.busy":"2023-09-24T15:30:58.607146Z","iopub.execute_input":"2023-09-24T15:30:58.607411Z","iopub.status.idle":"2023-09-24T15:30:58.615244Z","shell.execute_reply.started":"2023-09-24T15:30:58.607387Z","shell.execute_reply":"2023-09-24T15:30:58.614078Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"gene_list = list(gene_expr[gene_expr < q0].index)\nif len(gene_list) > 150:\n    gene_list = gene_list[:100]\ndownreg = enrichr_results(gene_list)\ndownreg.results.head(10)","metadata":{"execution":{"iopub.status.busy":"2023-09-24T15:30:58.616770Z","iopub.execute_input":"2023-09-24T15:30:58.617314Z","iopub.status.idle":"2023-09-24T15:31:01.363294Z","shell.execute_reply.started":"2023-09-24T15:30:58.617276Z","shell.execute_reply":"2023-09-24T15:31:01.362277Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"try:\n    dot_plot(downreg, cutoff=0.10)\nexcept ValueError:\n    print(\"No pathway match was significant for this gene list!\")","metadata":{"execution":{"iopub.status.busy":"2023-09-24T15:31:01.364398Z","iopub.execute_input":"2023-09-24T15:31:01.364718Z","iopub.status.idle":"2023-09-24T15:31:04.000140Z","shell.execute_reply.started":"2023-09-24T15:31:01.364682Z","shell.execute_reply":"2023-09-24T15:31:03.999162Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Exploring a Limited Gene Set + Compound Class Expression Space\n\nMy prior EDA on both [chemical structure and supplementary metadata](https://www.kaggle.com/code/vendekagonlabs/supplemental-compound-info/) and GSEA above are both in service of one goal: Can we limit the prediction problem to: certain categories of compounds impact certain pathways in predictable ways, and consider the rest as noise? (The rest is most likely _not_ noise, but this is a simplifying fiction for some exploratory work).\n\nTo do that, in the section below I read in the supplemental data added and cleaned from the notebook linked in the previous paragraph, and read in a well known and studied pathway (mTOR) from KEGG, and then use those to explore a small subset of drugs expected to impact mTOR genes.","metadata":{}},{"cell_type":"code","source":"supp_data_dir = Path('/kaggle/input/supplemental-compound-info')\ncompound_df = pd.read_csv(supp_data_dir / 'compounds.tsv', sep='\\t')\nprint(compound_df.columns)\ncompound_df.head(5)","metadata":{"execution":{"iopub.status.busy":"2023-09-24T15:31:04.001362Z","iopub.execute_input":"2023-09-24T15:31:04.001691Z","iopub.status.idle":"2023-09-24T15:31:04.031362Z","shell.execute_reply.started":"2023-09-24T15:31:04.001658Z","shell.execute_reply":"2023-09-24T15:31:04.030367Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We read in information on the mTOR signaling pathway from the Maayan Lab, whose web infrastructure backs a lot of the requests made by GSEA packages.","metadata":{}},{"cell_type":"code","source":"import requests\n\nmtor_url = \"https://maayanlab.cloud/Harmonizome/api/1.0/gene_set/mTOR+signaling+pathway/PID+Pathways\"\nresp = requests.get(mtor_url)\nassert resp.status_code == 200","metadata":{"execution":{"iopub.status.busy":"2023-09-24T15:31:04.032539Z","iopub.execute_input":"2023-09-24T15:31:04.032853Z","iopub.status.idle":"2023-09-24T15:31:04.330946Z","shell.execute_reply.started":"2023-09-24T15:31:04.032828Z","shell.execute_reply":"2023-09-24T15:31:04.330046Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The results are returned in a nested json data structure, but we just need gene symbol names, which we extract as a list using a comprehension.","metadata":{}},{"cell_type":"code","source":"obj = resp.json()\nmtor_gene_set = [o['gene']['symbol'] for o in obj['associations']]\nprint(f\"Found {len(mtor_gene_set)} genes associated with the mTOR pathway.\")\nmtor_gene_set[:15]","metadata":{"execution":{"iopub.status.busy":"2023-09-24T15:31:04.332223Z","iopub.execute_input":"2023-09-24T15:31:04.332552Z","iopub.status.idle":"2023-09-24T15:31:04.341992Z","shell.execute_reply.started":"2023-09-24T15:31:04.332525Z","shell.execute_reply":"2023-09-24T15:31:04.340909Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We then subset our `de_train` DataFrame to contain only mTOR gene expression data. Then we select a few metadata columns from the original `de_train` to join back in.","metadata":{}},{"cell_type":"code","source":"mtor_expr = de_train[[col for col in de_train.columns if col in mtor_gene_set]]\nprint(f\"Matched {len(mtor_expr.columns)} genes from the gene expression data with mTOR genes.\")\nmtor_expr.head(5)","metadata":{"execution":{"iopub.status.busy":"2023-09-24T15:31:04.343290Z","iopub.execute_input":"2023-09-24T15:31:04.343687Z","iopub.status.idle":"2023-09-24T15:31:04.397220Z","shell.execute_reply.started":"2023-09-24T15:31:04.343655Z","shell.execute_reply":"2023-09-24T15:31:04.396403Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"mtor_subset = de_train[['sm_name', 'cell_type']].join(mtor_expr)\nmtor_subset.head(5)","metadata":{"execution":{"iopub.status.busy":"2023-09-24T15:31:04.398480Z","iopub.execute_input":"2023-09-24T15:31:04.398821Z","iopub.status.idle":"2023-09-24T15:31:04.426473Z","shell.execute_reply.started":"2023-09-24T15:31:04.398798Z","shell.execute_reply":"2023-09-24T15:31:04.425485Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now that we've isolated mTOR specific gen expression, we want to define a set of compounds known to inhibit mTOR activity. Most of these will be kinase inhibitors of one form or another.","metadata":{}},{"cell_type":"code","source":"compound_df.dropna(subset=['moa', 'target'], inplace=True)\ninhibitor_df = compound_df[compound_df['moa'].str.contains('inhibitor')]\ninhibitor_df.head(5)","metadata":{"execution":{"iopub.status.busy":"2023-09-24T15:31:04.427456Z","iopub.execute_input":"2023-09-24T15:31:04.427743Z","iopub.status.idle":"2023-09-24T15:31:04.450453Z","shell.execute_reply.started":"2023-09-24T15:31:04.427719Z","shell.execute_reply":"2023-09-24T15:31:04.449555Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We limit our metadata to only contain target, `sm_name` (the chemical name that joins to our gene expression DataFrame), and the targets and mechanism of action of the drug. Note the duplicates remaining here mean that a compound has more than one mechanism of action or more than one drug. This won't impact out downstream analysis as we'll just use the compound list from here on out.","metadata":{}},{"cell_type":"code","source":"mtor_inhibitor_df = inhibitor_df[inhibitor_df['target'].isin(mtor_gene_set)]\nprint(f\"{len(mtor_inhibitor_df['cmap_name'].unique())} compounds are \" + \\\n      \"categorized as inhibitors and target mTOR pathway genes.\")\nmtor_inhibitor_df = mtor_inhibitor_df[['target', 'sm_name', 'moa']].drop_duplicates()\nmtor_inhibitor_df","metadata":{"execution":{"iopub.status.busy":"2023-09-24T15:31:04.453414Z","iopub.execute_input":"2023-09-24T15:31:04.453744Z","iopub.status.idle":"2023-09-24T15:31:04.469995Z","shell.execute_reply.started":"2023-09-24T15:31:04.453716Z","shell.execute_reply":"2023-09-24T15:31:04.469119Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Visualizing mTOR inhibiting compounds vs mTOR pathway gene expression with UMAP\n\nNow that we've isolated a set of compounds that are listed as inhibitors and target genes in the mTOR pathway, we will look to see if they cluster together in the mTOR pathway expression space after dimensionality reduction.","metadata":{}},{"cell_type":"code","source":"mtor_inhibitor_sm_names = list(mtor_inhibitor_df['sm_name'].unique())\n\ndef mtor_category(sm_name):\n    if sm_name in mtor_inhibitor_sm_names:\n        return \"mtor_inhibitor\"\n    else:\n        return \"other\"\n\nmtor_subset['mtor_category'] = mtor_subset['sm_name'].apply(mtor_category)\n\n# sanity check that all 10 'mtor_inhibitor' compounds joined successfully\nmtor_subset[['sm_name', 'mtor_category']].drop_duplicates()['mtor_category'].value_counts()","metadata":{"execution":{"iopub.status.busy":"2023-09-24T15:31:04.474920Z","iopub.execute_input":"2023-09-24T15:31:04.475440Z","iopub.status.idle":"2023-09-24T15:31:04.486802Z","shell.execute_reply.started":"2023-09-24T15:31:04.475416Z","shell.execute_reply":"2023-09-24T15:31:04.485748Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## UMAP Visualization\n\nNow that we've got the metadata of interest annotated in our limited geneX space, we visualized the UMAP embeddings. To keep things simple, we first do this by limiting the cell type to CD4+ T cells.","metadata":{}},{"cell_type":"code","source":"%%capture\n# ^ suppress numba deprecation warnings for UMAP\nfrom umap.umap_ import UMAP\nimport plotly.express as px","metadata":{"execution":{"iopub.status.busy":"2023-09-24T15:35:04.654466Z","iopub.execute_input":"2023-09-24T15:35:04.654921Z","iopub.status.idle":"2023-09-24T15:35:04.662333Z","shell.execute_reply.started":"2023-09-24T15:35:04.654889Z","shell.execute_reply":"2023-09-24T15:35:04.661344Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"mtor_subset.columns","metadata":{"execution":{"iopub.status.busy":"2023-09-24T15:31:04.495849Z","iopub.execute_input":"2023-09-24T15:31:04.496363Z","iopub.status.idle":"2023-09-24T15:31:04.507696Z","shell.execute_reply.started":"2023-09-24T15:31:04.496330Z","shell.execute_reply":"2023-09-24T15:31:04.506204Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cd4 = mtor_subset[mtor_subset['cell_type'] == 'T cells CD4+'].copy(deep=True)\ngX = cd4[cd4.columns[2:-1]].values.copy()\nump = UMAP(n_components=3)\ngX_proj = ump.fit_transform(gX)\ngX_proj.shape","metadata":{"execution":{"iopub.status.busy":"2023-09-24T15:31:04.509119Z","iopub.execute_input":"2023-09-24T15:31:04.509546Z","iopub.status.idle":"2023-09-24T15:31:07.203192Z","shell.execute_reply.started":"2023-09-24T15:31:04.509512Z","shell.execute_reply":"2023-09-24T15:31:07.202120Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cd4['umapX'] = gX_proj[:, 0]\ncd4['umapY'] = gX_proj[:, 1]\ncd4['umapZ'] = gX_proj[:, 2]\ncd4['size_hint'] = cd4['sm_name'].apply(lambda x: 8 if x in mtor_inhibitor_sm_names else 2)","metadata":{"execution":{"iopub.status.busy":"2023-09-24T15:31:07.204446Z","iopub.execute_input":"2023-09-24T15:31:07.204784Z","iopub.status.idle":"2023-09-24T15:31:07.213325Z","shell.execute_reply.started":"2023-09-24T15:31:07.204754Z","shell.execute_reply":"2023-09-24T15:31:07.212547Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"px.scatter_3d(\n    cd4,\n    x='umapX',\n    y='umapY',\n    z='umapZ',\n    color='mtor_category',\n    size='size_hint',\n    hover_data=['sm_name']\n)","metadata":{"execution":{"iopub.status.busy":"2023-09-24T15:31:07.214709Z","iopub.execute_input":"2023-09-24T15:31:07.215025Z","iopub.status.idle":"2023-09-24T15:31:07.301209Z","shell.execute_reply.started":"2023-09-24T15:31:07.214998Z","shell.execute_reply":"2023-09-24T15:31:07.300231Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"These drugs seem to be distributed across the UMAP projection space, rather than grouped. It's possible there might be subgroupings, but we would need to find some criteria to validate this and rule out the more conservative working assumption that these are basically randomly distributed.\n\nContinuing, we wrap our visualization into one function so we can repeat it simply for other cell types.","metadata":{}},{"cell_type":"code","source":"def mtor_umap(df, limit_cell_types=None):\n    \"\"\"Given a limited gene expression DataFrame `df` and an optional\n    list of `limit_cell_types` to restrict cell types used, preforms\n    dimensionality reduction w/UMAP and visualizes the resulting embeddings\n    in a 3d plot, with mtor inhibitors colored and larger.\"\"\"\n    if limit_cell_types is not None:\n        sub = df[df['cell_type'].isin(limit_cell_types)].copy(deep=True)\n    else:\n        sub = df.copy(deep=True)\n    gX = sub[sub.columns[2:-1]].values.copy()\n    ump = UMAP(n_components=3)\n    gX_proj = ump.fit_transform(gX)\n    gX_proj.shape\n    sub['umapX'] = gX_proj[:, 0]\n    sub['umapY'] = gX_proj[:, 1]\n    sub['umapZ'] = gX_proj[:, 2]\n    sub['size_hint'] = sub['sm_name'].apply(lambda x: 8 if x in mtor_inhibitor_sm_names else 2)\n    return px.scatter_3d(\n        sub,\n        x='umapX',\n        y='umapY',\n        z='umapZ',\n        color='mtor_category',\n        size='size_hint',\n        hover_data=['sm_name', 'cell_type']\n    )","metadata":{"execution":{"iopub.status.busy":"2023-09-24T15:31:07.302740Z","iopub.execute_input":"2023-09-24T15:31:07.303494Z","iopub.status.idle":"2023-09-24T15:31:07.310905Z","shell.execute_reply.started":"2023-09-24T15:31:07.303463Z","shell.execute_reply":"2023-09-24T15:31:07.310179Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"mtor_subset['cell_type'].unique()","metadata":{"execution":{"iopub.status.busy":"2023-09-24T15:31:07.311965Z","iopub.execute_input":"2023-09-24T15:31:07.312410Z","iopub.status.idle":"2023-09-24T15:31:07.328120Z","shell.execute_reply.started":"2023-09-24T15:31:07.312385Z","shell.execute_reply":"2023-09-24T15:31:07.327241Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Here we'll limit the gene expression space to all T cells, then only NK cells.","metadata":{}},{"cell_type":"code","source":"mtor_umap(mtor_subset,\n          limit_cell_types=['T cells CD4+',\n                            'T cells CD8+',\n                            'T regulatory cells'])","metadata":{"execution":{"iopub.status.busy":"2023-09-24T15:31:07.329138Z","iopub.execute_input":"2023-09-24T15:31:07.329423Z","iopub.status.idle":"2023-09-24T15:31:10.594205Z","shell.execute_reply.started":"2023-09-24T15:31:07.329399Z","shell.execute_reply":"2023-09-24T15:31:10.593246Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"mtor_umap(mtor_subset,\n          limit_cell_types=['NK cells'])","metadata":{"execution":{"iopub.status.busy":"2023-09-24T15:31:10.595670Z","iopub.execute_input":"2023-09-24T15:31:10.596041Z","iopub.status.idle":"2023-09-24T15:31:13.354440Z","shell.execute_reply.started":"2023-09-24T15:31:10.596006Z","shell.execute_reply":"2023-09-24T15:31:13.353420Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# More coming\n\nThe next section I'll be adding will use GSEA methods on the single cell data in subsets. I may put this into another notebook so that this one will continue to load fast, allowing faster exploration of simple enrichment with the main training dataset.\n\nIf you found this notebook helpful, please upvote and check out the other notebooks I've contributed to this contest:\n\n- [Single Cell EDA with Scanpy](https://www.kaggle.com/code/vendekagonlabs/op2-single-cell-eda-10x-multiome)\n- [Computational Chemistry with our Small Molecules in R](https://www.kaggle.com/code/vendekagonlabs/come-on-chemicals-r-version)\n- [Adding Supplemental Compound Info](https://www.kaggle.com/code/vendekagonlabs/supplemental-compound-info/)\n","metadata":{}}]}