{"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":"# Exploring gene overlap between Broad LINCS data and the competition dataset\n#### Part 1 of (hopefully) a series exploring the massive Broad LINCS dataset and explaining the basic concepts\n\nThe Broad LINCS dataset contains 12328 rows (genes) by 118050 columns (samples) for a total of 1,455,320,400 entries\n\nIn this notebook we will be examining the metadata about the genes, contained in the file `GSE70138_Broad_LINCS_gene_info_2017-03-06.txt`\n\nThe Broad LINCS dataset, including the metadata is available on the [GEO Website](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE70138). If there is any terminology I use that you are not familiar with, it is probably in the [Connectopedia Glossary](https://clue.io/connectopedia/glossary).","metadata":{}},{"cell_type":"code","source":"import pandas as pd\n# All metadata files are tab delimited text files\ngene_df = pd.read_csv(\"/kaggle/input/gse70138-broad-lincs-gene-info-2017-03-06/GSE70138_Broad_LINCS_gene_info_2017-03-06.txt\", sep=\"\\t\")\n# Count the number of metadata entries.\ngene_df.shape[0]","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-09-29T16:23:23.624735Z","iopub.execute_input":"2023-09-29T16:23:23.625124Z","iopub.status.idle":"2023-09-29T16:23:23.652960Z","shell.execute_reply.started":"2023-09-29T16:23:23.625097Z","shell.execute_reply":"2023-09-29T16:23:23.651892Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As expected, we have one metadata entry, for every row in the LINCS data.\nLet's see what metadata is available about genes:","metadata":{}},{"cell_type":"code","source":"gene_df.columns","metadata":{"execution":{"iopub.status.busy":"2023-09-29T16:23:23.654488Z","iopub.execute_input":"2023-09-29T16:23:23.654812Z","iopub.status.idle":"2023-09-29T16:23:23.660659Z","shell.execute_reply.started":"2023-09-29T16:23:23.654789Z","shell.execute_reply":"2023-09-29T16:23:23.659543Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The one's we are interested in are:\n#### `pr_gene_symbol`, the symbolic name of the gene.\nThese correspond to the column names in the training data.","metadata":{}},{"cell_type":"code","source":"gene_df[\"pr_gene_symbol\"][:10]","metadata":{"execution":{"iopub.status.busy":"2023-09-29T16:23:23.662059Z","iopub.execute_input":"2023-09-29T16:23:23.662383Z","iopub.status.idle":"2023-09-29T16:23:23.674053Z","shell.execute_reply.started":"2023-09-29T16:23:23.662360Z","shell.execute_reply":"2023-09-29T16:23:23.673200Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### `pr_is_lm`, whether the gene is a landmark gene.\n\nBecause the number of genes to measure is so high, scientists directly overserve a smaller set of genes and then use those observations to impute (predict) the data for the rest of the genes.","metadata":{}},{"cell_type":"code","source":"# Stored as an integer 1 or 0, instead of binary True or False.\ngene_df[[\"pr_gene_symbol\",\"pr_is_lm\"]][:10]","metadata":{"execution":{"iopub.status.busy":"2023-09-29T16:23:23.675801Z","iopub.execute_input":"2023-09-29T16:23:23.676409Z","iopub.status.idle":"2023-09-29T16:23:23.689199Z","shell.execute_reply.started":"2023-09-29T16:23:23.676386Z","shell.execute_reply":"2023-09-29T16:23:23.688537Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Further reading on on Landmark genes:\n* [\"What are landmark genes?\"](https://clue.io/connectopedia/what_are_landmark_genes)\n* [\"How were the 1000 landmarks identified?\"](https://clue.io/connectopedia/category/Analytical%20Methods)\n\nBy the way, some genes are easier to predict than others using the set of landmark genes. The easy set of genes is marked by `pr_is_bing` where BING stands for \"best inferred genes\".\n\n*Note*: Landmark genes are marked as BINGs, even though they are measured and not inferred.","metadata":{}},{"cell_type":"code","source":"gene_df[[\"pr_gene_symbol\",\"pr_is_lm\",\"pr_is_bing\"]][:10]","metadata":{"execution":{"iopub.status.busy":"2023-09-29T16:23:23.690046Z","iopub.execute_input":"2023-09-29T16:23:23.690914Z","iopub.status.idle":"2023-09-29T16:23:23.703248Z","shell.execute_reply.started":"2023-09-29T16:23:23.690880Z","shell.execute_reply":"2023-09-29T16:23:23.702675Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\"*Cool, but how does that help us in the competition?*\"\n\nGlad you asked! As it turns out, the competition dataset also contains measurements for the landmark genes! Let's take a look:\n\n*Note: The observations in the training data are transposed relative to the observations in the LINCS. In the training data, columns are genes, and rows are measurements.","metadata":{}},{"cell_type":"code","source":"# Load in training data.\ntrain_df = pd.read_parquet(\"/kaggle/input/open-problems-single-cell-perturbations/de_train.parquet\")\ntrain_genes = train_df.columns[5:]\ntrain_genes","metadata":{"execution":{"iopub.status.busy":"2023-09-29T16:23:23.704030Z","iopub.execute_input":"2023-09-29T16:23:23.704623Z","iopub.status.idle":"2023-09-29T16:23:24.782259Z","shell.execute_reply.started":"2023-09-29T16:23:23.704598Z","shell.execute_reply":"2023-09-29T16:23:24.781196Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's see which genes these datasets have in common:","metadata":{}},{"cell_type":"code","source":"# Find all genes (by index) in the metadata where the gene symbol is in the training data.\nshared_gene_idxs = gene_df[\"pr_gene_symbol\"].isin(train_genes)\n# Select just the relevant columns in the metadata for those indexes.\nshared_genes = gene_df[shared_gene_idxs][[\"pr_gene_symbol\",\"pr_is_lm\",\"pr_is_bing\"]]\nshared_genes","metadata":{"execution":{"iopub.status.busy":"2023-09-29T16:23:24.783550Z","iopub.execute_input":"2023-09-29T16:23:24.784229Z","iopub.status.idle":"2023-09-29T16:23:24.799537Z","shell.execute_reply.started":"2023-09-29T16:23:24.784200Z","shell.execute_reply":"2023-09-29T16:23:24.798547Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"*That's a lot of shared genes!*\n\nAs a quick check, let's see if every landmark gene in the LINCs data is also measured in the training data:","metadata":{}},{"cell_type":"code","source":"# Selected only the landmarks.\nlandmark_df = gene_df.loc[gene_df[\"pr_is_lm\"]]\n# Check that all landmark genes are in the training columns.\nlandmark_df[\"pr_gene_symbol\"].isin(train_genes).all()","metadata":{"execution":{"iopub.status.busy":"2023-09-29T16:23:24.801451Z","iopub.execute_input":"2023-09-29T16:23:24.802014Z","iopub.status.idle":"2023-09-29T16:23:24.816285Z","shell.execute_reply.started":"2023-09-29T16:23:24.801988Z","shell.execute_reply":"2023-09-29T16:23:24.815416Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"That's great!\n\n**If we can train a model to predict these landmark genes on the LINCS data, it shouldn't be too hard to transfer it to the competition data.**\n\n*Left as an exercise to reader :)*\n\nI wonder if all the BINGs are in the training data too?","metadata":{}},{"cell_type":"code","source":"# For every gene where 'pr_is_bing' is 1, check that the 'pr_gene_symbol' is in the train genes.\ngene_df.loc[gene_df[\"pr_is_bing\"]][\"pr_gene_symbol\"].isin(train_genes).all()","metadata":{"execution":{"iopub.status.busy":"2023-09-29T16:23:24.817582Z","iopub.execute_input":"2023-09-29T16:23:24.818245Z","iopub.status.idle":"2023-09-29T16:23:24.828940Z","shell.execute_reply.started":"2023-09-29T16:23:24.818211Z","shell.execute_reply":"2023-09-29T16:23:24.827876Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"That's a huge overlap. Let's visualize.","metadata":{}},{"cell_type":"code","source":"import seaborn as sns\nlincs_lincs_overlap = 1\nlincs_train_overlap = len(shared_genes)/len(gene_df)\ntrain_lincs_overlap = len(shared_genes)/len(train_genes)\ntrain_train_overlap = 1\noverlap_data = df=pd.DataFrame({\"NIH LINCS\":[lincs_lincs_overlap,lincs_train_overlap],\n                               'Training Data':[train_lincs_overlap,train_train_overlap]},index=['NIH LINCS','Training Data'])\n\nsns.heatmap(overlap_data,annot=True)","metadata":{"execution":{"iopub.status.busy":"2023-09-29T16:34:13.584428Z","iopub.execute_input":"2023-09-29T16:34:13.585104Z","iopub.status.idle":"2023-09-29T16:34:13.822889Z","shell.execute_reply.started":"2023-09-29T16:34:13.585074Z","shell.execute_reply":"2023-09-29T16:34:13.821908Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**This means that if you train on the shared genes in the NIH LINCS data, you've already covered half the targets in the competition data!**\n\nThat's all for this notebook. If you have any question, please feel free to leave a comment.\n\n#### Best of luck!","metadata":{}}]}